CC BY
131
Вестн. Сам. гос. техн. ун-та. Сер. Физ.-мат. науки. 2018. Т. 22, № 4. С. 669-701 ISSN: 2310-7081 (online), 1991-8615 (print) d http://doi.org/10.14498/vsgtu1643
Математическое моделирование, численные методы и комплексы программ
УДК 519.654
Численный метод нелинейного оценивания на основе разностных уравнений
В. Е. Зотеев
Самарский государственный технический университет, Россия, 443100, Самара, ул. Молодогвардейская, 244.
Аннотация
Рассматривается новый численный метод оценки параметров нелинейных математических моделей, в основе которого лежат разностные уравнения, описывающие результаты наблюдений. Алгоритм численного метода содержит следующие шаги:
- построение линейно-параметрической дискретной модели исследуемого процесса в форме разностных уравнений, коэффициенты которых известным образом связаны с параметрами нелинейной математической модели;
- формирование на основе разностных уравнений обобщенной регрессионной модели;
- вычисление оценки начального приближения и уточнения среднеквадратичных оценок коэффициентов обобщенной регрессионной модели на основе итерационной процедуры;
- вычисление оценок параметров нелинейной математической модели на основе среднеквадратичных оценок коэффициентов разностных уравнений;
- оценка погрешности результатов вычислений на основе методов статистической обработки данных эксперимента.
Предлагаются различные подходы к построению систем разностных уравнений для математических моделей в форме нелинейных функциональных зависимостей. Получены соотношения, лежащие в основе итерационного процесса уточнения коэффициентов обобщенной регрессионной модели, построенной на основе разностных уравнений. Описана процедура оценки погрешности результатов вычислений параметров нелинейных функциональных зависимостей, известным образом связанных с коэффициентами системы разностных уравнений. Применение
Научная статья
3 ©® Контент публикуется на условиях лицензии Creative Commons Attribution 4.0 International (https://creativecommons.org/licenses/by/4.0/deed.ru) Образец для цитирования
Зотеев В. Е. Численный метод нелинейного оценивания на основе разностных уравнений // Вестн. Сам. гос. техн. ун-та. Сер. Физ.-мат. науки, 2018. Т. 22, № 4. С. 669-701. doi: 10.14498/vsgtu1643. Сведения об авторе
Владимир Евгеньевич Зотеев А http://orcid.org/0000-0001-7114-4894 доктор технических наук, доцент; профессор каф. прикладной математики и информатики; e-mail: [email protected]
численного метода на основе разностных уравнений проиллюстрировано на примерах оценки параметров математической модели линейного осциллятора с затуханием, модели свободных колебаний диссипативной механической системы с турбулентным трением, а также параметров логистического тренда, описываемого функцией Верхулста (Перла-Рида).
Ключевые слова: математическая модель, нелинейный регрессионный анализ, система разностных уравнений, обобщенная регрессионная модель, среднеквадратическое оценивание, статистическая обработка результатов эксперимента.
Получение: 14 сентября 2018 г. / Исправление: 5 ноября 2018 г. / Принятие: 12 ноября 2018 г. / Публикация онлайн: 18 декабря 2018 г.
Введение. Одной из основных проблем при построении математических моделей по результатам наблюдений, полученных в ходе эксперимента или натурных испытаний, является нелинейность математической модели относительно ее параметров. В тех случаях, когда форма математической модели однозначно не обусловлена априорной информацией, полученной на основе физических, химических или иных законов, описанных, например, в виде уравнений математической физики, широко используются линейные регрессионные модели, оценивание и статистический анализ параметров которых, как правило, не вызывает каких-либо затруднений [1—3].
Однако достаточно часто математические модели в форме функциональных зависимостей, описывающих исследуемый процесс или взаимосвязь между характеристиками объекта исследования, строятся на основе решения систем интегро-дифференциальных уравнений, нелинейных по своей природе. Причем нередко линеаризация полученных соотношений настолько упрощает математическую модель, что применение ее в задаче параметрической идентификации становится не только нецелесообразным, но и практически бессмысленным. Таким образом, проблема нелинейного оценивания — вычисления параметров нелинейных математических моделей по результатам наблюдений — всегда была и остается важнейшей проблемой при математическом моделировании.
Пусть математическая модель наблюдаемого процесса описывается нелинейной функциональной зависимостью вида у = f (1,а\,а2,... ,ап), где аг, а2, ..., ап — параметры модели; п — число параметров; £ е [0, те) — аргумент функциональной зависимости, в частности, имеющий размерность времени. При равномерной дискретизации с периодом т: = тк, к = 0,1, 2,..., отсюда получаем дискретную математическую модель вида у = f (тк, аг,а2,..., ап), к = 0,1, 2, 3 ..., нелинейную по параметрам аг, а2,..., ап.
Пусть у к — результаты наблюдений (например, данные эксперимента), к = 0,1,2,3,...,Ы — 1; N — объем выборки результатов наблюдений. Тогда математическую модель, описывающую результаты наблюдений, можно представить в виде
Ук = / (тк, а\,а2,..., ап) + ек, к = 0,1,2,3,...,Ж — 1, (1)
где — разброс результатов наблюдений относительно модели ук при заданных значениях ее параметров (остатки). В векторной форме модель (1) имеет
следующий вид:
У = /( а) + £,
где у е , у = / (а) е Мга, ее , а = (аъа2,..., ап)т е Мга.
Уравнение (1) описывает регрессионную модель, нелинейную по своим параметрам а\, а2,..., ап. В качестве критерия близости модели (1) результатам наблюдений у^ в регрессионном анализе принимается норма разности в евклидовом пространстве N-мерных векторов [2,3]:
\
N-1
Yl(ук - Ук )2
к=0
где у, у G Е^ со скалярным произведением
\\у - У\\ = (У - У^ - У) = (У - у)Т(у - у).
С учетом выбранного критерия близости параметры нелинейной регрессионной модели (1) находятся из минимизации нелинейной функции
N -1
Q(ai,a2, ...,ап) = \\е\\2 = \\у - у\\2 = ^ [ук - f(r к, äi, а2,..., än)] ^ min.
к=о
В такой постановке эта задача относится к задаче нелинейного оценивания и может быть решена методами нелинейной регрессии, описанными в [1-6].
1. Классические методы нелинейного оценивания [2—9]. Известные традиционно применяемые в практике построения математических моделей методы нелинейного оценивания можно разбить на три основные группы. К первой группе относятся численные методы поиска минимума функции нескольких переменных, среди которых наибольшее применение нашли градиентный метод, метод наискорейшего спуска и метод сопряженных градиентов [3,10]. В основе этих методов лежит вычисление градиента функции Q(ä):
gradQ(ä) = -öQ = —2WT[y - f(ä)],
где Wт — матрица Якоби размера [Nxn], элементы которой = —
значения производных функции f(ä) по параметрам äj, j = 1, n, вычисленные в точках tк, к = 0,1, 2,..., N — 1.
Метод скорейшего спуска отличается от градиентного метода оптимальным выбором шага на каждой итерации. В свою очередь метод сопряженных направлений отличается от метода наискорейшего спуска только выбором направления уменьшения функции на каждом шаге — вектор
■p(i) = gradQ(a(i)) + a(i)p(i-1)
вместо вектора градиента gradQ(a(t)). Алгоритм метода сопряженных направлений может быть описан следующей системой соотношений:
а(°) =0, р(°) = цгаё((а(°)), а« = 1 2,
^гаё^а^ )|2
р(г) = ёга<1((а('1))+а(-'1)р(-'1-1), ^ = а^тт(а« — р • Р(г)), а (*+г) = а(г) — ¿V0, ¿ = 0,1,2, 3,...
При а(г) = 0 получаем формулы, описывающие алгоритм метода наискорейшего спуска. Если при этом шаг на каждой итерации задается некоторым произвольным образом, то приходим к градиентному методу. Основным недостатком методов непосредственного поиска минимума функции является их медленная сходимость и «зацикливание» в окрестности точки минимума.
Ко второй группе относятся методы, ориентированные непосредственно на решение нелинейной системы нормальных уравнений:
9( о) = ^ = —2^ (а)т[у — Да)]=0.
В основе методов этой группы лежит метод Ньютона [2,3]. Алгоритм вычислений по методу Ньютона описывается итерационной формулой
а(ш) = а« + 2Р(а(^)-1Ш(аЩт [у — /(а®)],
где
N -1
Р(а) = 2Ж(а)тЖ(а) — 2 ^ (ук — Д(а))Нк(а)
к=0
— матрица Якоби для вектор-функции д(а) е Мга размера [пхп], Нк(а) — матрицы Гессе, размера [пхп], элементы которых ., г,] = 1,п, вычисляются в точке í к ,к = 0,1, 2,..., N — 1. * '
В нелинейной регрессии обычно используют метод Ньютона—Гаусса, который отличается от метода Ньютона тем, что при нахождении матрицы Яко-би в формуле Р(а) пренебрегают слагаемыми, содержащими матрицы Гессе Нк(а). Итерационная процедура уточнения оценок коэффициентов нелинейной регрессионной модели методом Ньютона—Гаусса описывается соотношением
а(т) = + (а,(^)тШ(а«)]- V(а,(^)т [у — /(а({))].
Модификациями метода Ньютона—Гаусса являются метод Хартли и метод Левенберга—Марквардта [3,8,9]. В методе Хартли для ускорения сходимости длина шага делается переменной за счет коэффициента > 0. Метод Левенберга—Марквардта целесообразно использовать при плохо обусловленной матрице Р(а) = 2Ш(а)тЖ(а). Этот метод имеет много общего с методом регуляризации (ridge-оценок), применяемом при решении плохо обусловленных систем линейных алгебраических уравнений. Алгоритм метода строится на основе метода Ньютона—Гаусса и описывается рекуррентной формулой
а (*+1) = а« + [Ж (а« )тЖ (а«) + р« Да«)] (а«)т [у — Да«)],
где ^ > 0 — малый параметр, = Е в методе Левенберга (Е — еди-
ничная матрица) и = Е ■ (а^^Ш(а(г))] в методе Марквардта.
Стратегия выбора параметра ^ описана в [3].
К третьей группе классических методов нелинейного оценивания можно отнести методы, в основе которых лежит линеаризация нелинейной математической модели [3]. Различают параметрическую линеаризацию нелинейной функциональной зависимости посредством разложения в ряд Тейлора и линеаризацию посредством некоторого подходящего преобразования нелинейной функции регрессоров (например, с помощью логарифмирования). В первом случае задача сводится к линейной регрессии, причем алгоритм вычислений описывается той же итерационной формулой, что и в методе Ньютона—Гаусса.
Завершая анализ известных методов нелинейного оценивания, можно сделать следующие выводы.
Во-первых, эти методы, как правило, требуют вычисления производных сложной нелинейной функции до второго порядка включительно.
Во-вторых, в методе Ньютона—Гаусса, а также при параметрической линеаризации на основе разложения в ряд Тейлора используется аппроксимация нелинейной регрессионной модели линейной, что приводит к методической погрешности в результатах вычислений и, как следствие, смещению оценок параметров модели.
В-третьих, основной проблемой при применении известных методов является проблема выбора начального приближения, обеспечивающего сходимость итерационной процедуры к точке минимума. Этот основной недостаток существенно ограничивает область применения методов нелинейной регрессии при построении математических моделей по результатам эксперимента.
Предлагаемый новый численный метод нелинейного оценивания в формате описанной выше систематизации может быть отнесен к третьей группе методов, так как в нем оценивание параметров нелинейной регрессии сводится к вычислению коэффициентов линейно параметрической дискретной модели (ЛПДМ) в форме разностных уравнений. Однако в отличие от параметрической линеаризации на основе разложения в ряд Тейлора, ЛПДМ в форме разностных уравнений эквивалентна дискретной нелинейной модели, и методическая погрешность, обусловливающая смещение оценок параметров модели, отсутствует. Вторым существенным достоинством предлагаемого метода является нахождение оценки начального приближения вектора коэффициентов ЛПДМ непосредственно в процессе реализации алгоритма вычислений; при этом находить оценку начального приближения вектора параметров нелинейной модели не нужно.
2. Алгоритм численного метода нелинейного оценивания на основе разностных уравнений. В основе численного метода нелинейного оценивания на основе разностных уравнений лежит переход от нелинейной по параметрам математической модели к линейно параметрической дискретной модели, коэффициенты которой известным образом связаны с искомыми параметрами. Построенная на основе ЛПДМ обобщенная регрессионная модель, описывающая результаты наблюдений, позволяет свести задачу нелинейного регрессии к вычислению среднеквадратичных оценок коэффициентов ЛПДМ, т.е. к задаче линейного прикладного регрессионного анализа, решение которой не вызывает особых затруднений.
Численный метод нелинейного оценивания на основе разностных уравнений включает следующие основные этапы [8]:
- построение математической модели исследуемого процесса в форме линейных разностных уравнений, коэффициенты которых известным образом связаны с параметрами нелинейной математической модели;
- формирование разностных уравнений, описывающих результаты наблюдений;
- формирование на основе разностных уравнений элементов обобщенной регрессионной модели;
- вычисление оценки начального приближения коэффициентов обобщенной регрессионной модели;
- итерационная процедура уточнения среднеквадратичных оценок коэффициентов обобщенной регрессионной модели;
- вычисление оценок параметров нелинейной математической модели на основе среднеквадратичных оценок коэффициентов разностных уравнений;
- оценка погрешности результатов вычислений на основе методов статистической обработки данных эксперимента.
Рассмотрим поэтапную реализацию описанного выше алгоритма численного метода на основе разностных уравнений.
3. Построение математической модели исследуемого процесса в форме разностных уравнений. Под математической моделью исследуемого процесса в форме разностных уравнений будем понимать линейную комбинацию функций ({/к, ЪаЪу^г,..., Ук-г), .] = 1, п, связывающих (г + 1) последовательных значений ук, ук-1, ..., Ук-г дискретной нелинейной модели Ук = у(тк, а,1, а,2,..., ап) с коэффициентами Л1, Л2, ..., Хп:
п
Ьк+1 (Ук, Ук-Ъ ..., Ук-Г) = (Ук, Ук-Ъ ..., Ук-Г), (2)
3=1
где к = 0,1, 2,...,Ы — 1. При к < г будем считать, что переменные Ук-Г в разностном уравнении (2) отсутствуют.
При известных соотношениях
Л1 = Л1(01,02,..., ап),
Л2 = Л2(Q'l,й2,... ,an), (3)
Лп = Лп(а1, а,2,..., ап),
связывающих параметры а = (а1, а,2,..., оп)т нелинейной регрессионной модели (1) с коэффициентами Л = (Л1, Л2,..., Лп)т разностного уравнения (2), задача нелинейного оценивания сводится к задаче линейной алгебры — решению нормальной системы линейных алгебраических уравнений относительно коэффициентов Л^, ] = 1,п.
Рассмотрим три основных подхода к построению разностных уравнений.
3.1. Построение системы разностных уравнений на основе г-преобразования. Этот подход к построению разностных уравнений целесообразно использовать в тех случаях, когда дискретную нелинейную зависимость можно представить в виде линейной комбинации элементов последова-
тельностей, для которых известны и приведены в соответствующей литературе формулы ^-преобразования [12]. При этом применение первой и второй теорем смещения, правил дифференцирования изображений, свертки последовательностей и т.п. позволяет упростить процесс формирования разностных уравнений.
Рассмотрим применение данной методики на конкретном примере построения системы разностных уравнений для нелинейной математической модели линейного осциллятора с затуханием, последовательность значений которой описывается дискретной функцией вида [11]
Ук = ао ехр(—атк) еов(штк + фо), к = 0,1, 2,..., (4)
содержащей четыре параметра: динамические характеристики а и ш и параметры ао и фо, связанные с начальными условиями.
С учетом известных формул ^-преобразования [9]:
и
имеем
, и . , -г а вш т
х{а^ вттк} =
х2 — 2а х еов т + а2 2
к х — х а еов т
х[а еовтк} =
г2 — 2а х еовт + а2
х{ у к } = х{ао ехр(—атк) еов(штк + фо)} = а0 еовфох {ехр(—атк) еовшт к}—
, , , . ,, х2 еовфо — хехр(—ат)еовшт еовфо
— ао втфох{ехр(—ат) вшштк} = ао^—--------—-—-■
,г2 — 2 ехр(—ат) х еов шт + ехр(—2ат)
аох ехр(—ат) вт шт вт фо
х2 — 2х ехр(—ат) еов шт + ехр(—2ат)
х2 еов фо — х ехр(—ат) еов шт еов(шт — фо)
ао
х2 — 2хехр(—ат) еовшт + ехр(—2ат)
Вводя обозначения А1 = 2ехр(—ат)еовшт и Х2 = — ехр(—2ат), из последнего получаем
(1 — А1х-1 — \2х-2)х{ук} = ао еовфо — х-1ао ехр(—ат) еов(шт — фо). Отсюда с учетом первой теоремы смещения [12]
г-г ■ *{ Ук } = Ук-т, г = 0,1,2,... и условии, что при к — г < 0 принимается ук-г = 0, имеем
Ук — А1 Ук-1 — \2Ук-2 = ао еовфо5к — ао ехр(—ат)еов(шт — фо)5к-1, 1, при к = 0,
где 5к = , _ , , _
0, при к = 0.
При к ^ 2 получаем разностное уравнение вида
ук = А1 Ук-1 +А2 Ук-2, к = 2, Ъ,...^ — 1.
К этим ( N — 2) уравнениям необходимо добавить еще два уравнения, соответствующих к = 0 и к = 1. Полагая к = 0, получаем у0 = а0 cosф0. При к = 1 соответственно получаем
У1 — Л1 уо = —ао ехр(—ат) еов(шт — фо)
или
ух = 2а0 exp(-ат) cos шт cos ф0 — а0 exp(-ат) cos шт cos ф0 —
— а0 exp(—ат) sin шт sin ф0 = а0 exp(—ат) cos(wt + ф0).
Вводя дополнительно два коэффициента А4 = а0 exp(—ат) cos(cc>T + ф0) и Аз = а0 cosф0, получаем линейно-параметрическую дискретную модель линейного осциллятора с затуханием в форме системы разностных уравнений:
( Уо = ^
Ух = А4, (5)
I Ук = Ахук-х + А2Ук-2, к = 2,3,4,...,
в которой коэффициенты Aj связаны с параметрами нелинейной зависимости (4) соотношениями
Ах = 2exp(—ат) cos шт, А2 = exp(—2ат),
А3 = а0 cos ф0, А4 = а0 exp(—ат)cos(í^T + ф0).
(6)
3.2. Построение системы разностных уравнений с использованием элементов матричной алгебры. Вначале рассмотрим применение данной методики к построению системы разностных уравнений на примере дискретной нелинейной модели, описывающей ординаты свободных колебаний диссипативной системы с турбулентным трением [11,13]:
Ук = ! ^ ^(штк + фо). (7)
Из соотношения (1 + атк)ук = а0 cos(c^Tk + ф0) получаем
[1 + ат (к — 1)]ук-х = а0 cos(^Tk + ф0) cos шт + а0 sin(штк + ф0) sin шт, [1 + ат (к — 2)]ук-2 = а0 cos(í^Tk + ф0) cos 2шт + а0 sin(штк + ф0) sin 2шт.
Обозначим Ак = а0 cos(штк + ф0), Вк = а0 sin(штк + ф0). Тогда получаем систему двух алгебраических уравнений относительно переменных Ак и Вк:
J Ак cosшт + Вк sinшт = [1 + ат (к — 1)}ук-х, 1 Ак cos 2шт + Вк sin 2шт = [1 + ат (к — 2)]ук-2,
из решения которой находим
Ак = 2 cos шт [1 + ат (к — 1)]ук-х — [1 + ат (к — 2)}ук-2.
Так как Ак = (1 + атк)ук, имеем
(1 + атк)ук = 2 cos шт [1 + ат (к — 1)]ук-х — [1 + ат (к — 2)]ук-2.
Отсюда, обозначая А1 = 2coswt и А2 = ат, можно получить разностное уравнение вида
Ук + Ук-2 = А1 Ук-1 - Ук - - 1)Ук-1 + (к - 2)ук-2], к = 2,3,....
Так как число параметров в нелинейной модели равно четырем: а = = (а,ш,а0,ф0)т, а число коэффициентов в системе разностных уравнений — двум, добавим в систему еще два уравнения с двумя коэффициентами Аз и A4, с учетом которых можно найти параметры ао и фо:
Уо = A3, yi = A4,
где
А3 = а0 еовф0, А4 =-cos(wt + ф0).
1 + ат
Таким образом, линейно-параметрическая модель, описывающая ординаты колебаний диссипативной системы с турбулентным трением, в форме системы разностных уравнений имеет вид
Уо = Аз,
Vi = А4, /g\
Ук + У к—2 = А1 Ук-1 - А2[кук - А1(к - 1)ук-1+ у J
+(к - 2)ук—2], к = 2, 3,...,
в которой коэффициенты связаны с параметрами нелинейной зависимости (7) соотношениями
А1 = 2coswT, А2 = ат, А3 = а0 соаф0, А4 = ——— cos(wt + ф0). (9)
1 + ат
Решим задачу построения системы разностных уравнений с использованием элементов матричной алгебры в общем виде для нелинейных математических моделей, дискретный аналог которых удовлетворяет следующему условию (достаточное условие):
m
дг[ук—г, т(к - г), а] fj (тк,а)фз (г, а), г = 0,1, 2,..., к^ г, (10)
3 = 1
где функции фз(г, а) не зависят от тк и ук—г; а = (а1, а2,..., ап)т — вектор неизвестных коэффициентов; m < п.
Полагая в (10) г = 1, 2, 3,... ,m, получаем систему линейных алгебраических уравнений относительно переменных fj(тк,а) (j = 1,m):
( ф1(1, а) fi(T к, а) + ф2(1, а) f2(тк, а)) + ■■■
+фт(1,а) fm(rк, а) = д1[ук—1, т(к - 1),а],
ф1(2, а) fi(T к, а) + ф2(2, а) f2(тк, а)) + ■■■
+фт(2,а)/т(тк,а) = д2[Ук—2, т(к - 2), а], (ц)
ф1(т, а)}\(тк, а) + ф2(т, а)¡2(тк, а)) +----
+фт(т, а) fm(rк, а) = дт[Ук—т, т(к - т), а];
Введем в рассмотрение векторы
/(тк, а) = [ Л(гк, а), ¡2(тк, а),..., /т(тк, а)]т,
(г, а) = [ф1 (г, а),ф2(г ,а),..., фт( г, а)\
т
V = (ук-1, Ук-2, . . ., Ук-т)т,
д(У, тк, а)= [§1 (ук-Ъ т(к—1), а), д2 (Ук-2, т(к—2),а),..., дт( ук-т, т(к—т), а)]
т
и матрицу системы линейных уравнений
Ф(а) = ф(1, а)
ф(2,а)
ф(т, а)
ф1(1,а) ф2(1,а) ф1(2,а) М2,а)
т
фт(1,а)
фт(2,а)
ф1(т,а) ф2(т,а) ... фт(т,а).
Тогда при условии ёе!Ф( а) = 0 решение системы линейных уравнений (11) можно представить в виде
¡(тк,а) = Ф(а)-1д (у,тк,а).
С учетом этой формулы получаем разностное уравнение вида
9о(Ук ,тк,а) = ф(0,а) Ф(а) ■ д(ук-1, Ук-2, ••• , Ук-т,тк,а), к ^ т. (12) Обозначая через Л = Ф(а)-1ф(0,а) вектор размерности [ т х 1], из ( ) имеем
9о(Ук ,г к, а) = Лт ■ д(ук-1, Ук-2,..., Ук-т,гк,а) =
т-> т
^ (ук-у, т(к — j),a), к ^ 3 = 1
т.
К этим уравнениям, число которых должно быть не менее т, содержащим т коэффициентов Л^, ] = 1, т, известным образом связанных с параметрами а = (а1,а2,..., ап)т:
Л = Ф(а)-1ф(0,а), можно добавить еще (п — т) уравнений
Уо = Лт+1, У1 = Лm+2,
Уп-т-1 — Лг.
содержащих (п — т) коэффициентов, также известным образом связанных с параметрами а = (а1, а2,..., ап)т:
Лт+1 = у(0, а,1,а,2,..., ап), Лт+2 = у(т,0,1,0,2,.. .,ап), Лт+з, = у(2т, 01, 02,..., ап),
, Лп = у[(п — т — 1) Т, а,1,а,2,..., ап],
где у(тк, а,1,а,2,..., оп) — известная нелинейная функциональная зависимость.
Для рассмотренного выше примера модели свободных колебаний систем с турбулентным трением (7) имеем
9г [ук-г, т(к - г), о] = [1 + ат(к - г)]ук-г, ¡\(тк, о) = о0 сов(штк + ф0), /2(тк, о) = о0 в'т(штк + ф0),
ф\(г,о) = СОИШТГ, ф2(г,о) = 8ШШТГ, ф(г, о) = (С008ШТГ, вШШТг)Т,
ф(0,о) = (1, 0)Т, ф(1,о) = (соишт, вшшт)Т, ф(2, о) = (сои2шт, вт2шт)Т.
Здесь т = 2, о = (о\, о2,..., оп)Т — вектор неизвестных параметров. Матрица Ф в этом случае имеет вид
Ф
cos шт sin шт cos 2шт sin 2шт
Определитель матрицы det Ф = sin шт = 0 при 0 < т < Т/2, где Т = 2ж/ш — период колебаний. Тогда обратную матрицу Ф-1 можно представить в виде
Ф
-1
1
sin ш
sin 2шт — sinшт cos 2 ш cos ш
2 cos ш 1
ctgшт
cos 2шт
При этом
ф(0,а)тФ(а)-1 = (1, 0) ■
2 cos ш — 1
— С0П2шт ^шт
= (2cosшт, — 1),
9o(ilk, тк, а) = (а + атк)ук,
д(ук-1,ук-2,тк, а) = ([1 + ат(к — 1)}ук-1, [1 + ат(к — 2)]ук-2) В соответствии с полученной выше формулой (12) имеем
т
до(ук, тк, а) = ф(0, а)тФ(а) 1 ■ д(ук-1, ук-2,тк, а)
или
(1 + атк)ук = (2cosшт, —1)
[1 + ат (к — 1)\ук-1 [1 + ат (к — 2)}ук-2
= 2 cos шт [1 + ат (к - ук-1 - [1 + ат (к - 2)] ук-2.
Вводя обозначения А1 = 2 cos шт и Х2 = ат, после простых преобразований получаем разностное уравнение
Ук - У к—2 = Аук-1 - МкУк - Ai(k - 1)ук-1 + (к - 2)ук-2], к = 2,3,...,
совпадающее с полученным выше результатом (8).
В частном случае, когда дг(Ук—г,тк,а) = ук—г, построенная система разностных уравнений принимает простой вид линейной комбинации отсчетов Ук—j, j = 1,ш, которую можно представить в форме скалярного произведения:
ук = f(Tк, а)т • ф(0, а) = ф(0, а)тФ-1 • Y
или
Ук = Ат - У = ^ AjУк-j, к ^т. з=1
Например, для рассмотренной выше дискретной модели линейного осциллятора с затуханием (6) имеем
Ук-Г = а0 ехр(—ат к) cos(шт к + ф0) ■ ехр(ат г) cos шт г+
+а0 exp(-атк) sin(wTк + фо) ■ ехр(атг) sin штг, f1(тк,а) = а0 ехр(-атк) cos(штk + ^0), ф1(г, а) = ехр(атг) cos штг, f2 (тк,а) = а0 ехр(-атк) sin(штк + ^0), ф2 (г, а) = ехр(ат г) sin шт г; ф1(0, а) = 1, ф1(1,а) = ехр(ат) cos шт, ф1(2,а) = ехр(2ат) cos 2шт, ф2(0,а) =0, ф2(1,а) = ехр(ат) sin шт, ф2(2,а) = ехр(2ат) sin 2шт.
Здесь т = 2, а = (а1, а2,..., ап)Т — вектор неизвестных параметров. Матрица Ф в этом случае имеет вид
Ф
ехр(ат) cos шт ехр(ат) sin шт ехр(2ат) cos 2шт ехр(2ат) sin 2 шт
= ехр( ат)
cos ш sin ш
ехр(ат) cos 2шт ехр( ат) sin 2шт
Определитель матрицы det Ф = exp(3 ат) sin шт = 0 при 0 < т < Т/2, где Т = 2тт/ш — период колебаний.
Тогда обратную матрицу Ф-1 можно представить в виде
Ф
1
1
ехр(3ат) sin шт
ехр(2 ат) sin 2шт — ехр( ат) sin шт — ехр(2ат) cos 2шт ехр( ат) cos шт
2ехр(—ат) cos шт — ехр(—2ат)
cos 2 ш cos ш — ехр(—ат)—- ехр(—2ат) —
sin ш
В соответствии с полученной выше формулой имеем
sin ш
Ук = (1, 0)
' 2ехр(—ат) cos шт — ехр(—2ат)
cos 2 ш cos ш
— ехр(—ат)- ехр(—2 ат)-
sin ш sin ш
= [2 ехр(—ат) cos шт, — ехр(—2ат)]
г -|
Ук-1
Ук-2_
Ук-1 Ук-2
= 2 exp(-ar) cos шт ■ ук-\ - exp(-2ar) ■ ук-2.
Вводя обозначения Х\ = 2exp(-ат) cos шт и Х2 = — exp(-2ат), получаем разностное уравнение второго порядка ук = Ук-\ + Х2 Ук-2, к = 2, 3,..., совпадающее с полученным выше результатом (5).
3.3. Построение системы разностных уравнений на основе формирования дифференциального уравнения и формул численного дифференцирования. В основе этого подхода лежит нахождение производных (обычно не более второго порядка включительно), формирование на их основе дифференциального уравнения с последующей заменой производных разностными отношениями в соответствии с известными формулами численного дифференцирования [10,14].
Рассмотрим применение данного подхода к формированию системы разностных уравнений на примере логистической функции Верхулста (Перла-Рида) [15-17]
№ =-^-г.
1 + а\ exp(-od)
Полагая tk = t0 + тк, к = 0,1, 2,..., где t0 — момент времени первого отсчета в результатах наблюдений, получаем дискретную трехпараметриче-скую нелинейную функцию вида
Ук =--г, к = 0,1, 2,..., (13)
1 + а\ exp(-octk)
для которой построим систему линейных разностных уравнений.
Дифференцируя обе части равенства y(t)+a\ exp(-at)y(t) = а0, получаем
y'(t) — аа\ exp(—at)y(t) + а\ exp(-od)y'(t) = 0. у'(t)[1 + а\ exp(—at)] = аа\ exp(—at)y'(t)
Отсюда
или
у (t)ао = оа^ exp(—at)y(t).
№
С учетом равенства
а\ ехр(—аЪ) = -т0- — 1
получаем дифференциальное уравнение вида
у'(1) = ау(1) — - у2(1). (14)
ао
К этому уравнению добавим еще одно уравнение, описывающее начальное условие и содержащее параметр а\:
У(*о) = -^-гг = Уо.
1 + а\ ехр(—а0)
Теперь, используя известные формулы численного дифференцирования (например, второго порядка точности [14]), дифференциальное уравнение (14) аппроксимируем системой разностных уравнений вида
—3уо + 4у1 — у2 „ а „2
---= аУо — — У1,
2т ао
Ук — Ук—2~ а „2
-—-= аУк-1 — — у2—, к = 0,1 2,...
¿т ао
Отсюда, обозначая
\ о \ 2та Х1 \ а° пк\ Хг = 2та, Х2 =--=--, Аз = —--, (15)
а0 а0 1 + а1 ехр(—а10)
с учетом начальных условий получаем систему разностных уравнений
Уо = Аз,
-3 Уо + 4 у 1 -У2 = Хгуо + Х2у0, (16)
Ук - Ук-2 = Хгук-1 + Х2 у1-г, к = 2,3,...,
в которой коэффициенты разностных уравнений связаны с параметрами функции Верхулста соотношениями (15).
Результаты численно-аналитических исследований показали, что методическая погрешность, связанная с аппроксимацией производной в (14) формулами численного дифференцирования, несущественна и находится в пределах разброса данных наблюдений относительно модели.
4. Построение разностных уравнений, описывающих результаты наблюдений, и формирование на их основе обобщенной регрессионной модели. Пусть данные наблюдений ук отличаются от результатов вычисления по модели у к на случайную величину естественного разброса £к:
Ук = Ук + £к, к = 0,1, 2,...,Ы - 1,
где N — объем выборки результатов наблюдений. С учетом этого равенства система разностных уравнений (2) принимает вид
Ьк+1 (Ук — £к, Ук-1 — £к-1, Ук-г — £к-г) =
п
У^ !'к+1^(Ук — £к, У к-1 — £к-1, Ук-г — £к-г), (17) 3=1
где к = 0,1, 2,...^ — 1, ук-г = 0 при к < г.
Используя разложение в общем случае нелинейных относительно Ук-] функций в (17) по степеням £к-], ограничиваясь при этом линейной составляющей, получаем
Ьк+1(Ук — £к, Ук-1 — £к-1,... , Ук-г — £к-г) ~
, , ч Ък+1(У к, Ук-Ъ..^ Ук-г) ~ Ьк+1(Ук, Ук-l,..., Ук-г) - -д-
1=0 °Ук-*
Iк+1^(Ук — £к, Ук-1 — £к-1,. . . , Ук-г — £к-г) ~
^ д!к+1,](Ук
д Ук-
, , ч ¡к+1,1(У к, Ук-Ъ..^ Ук-г) 1к+1,]{ук,Ук-l,...,Ук-г) - -—:-£к-
=о
С учетом этих соотношений аппроксимация системы разностных уравнений (17) может быть представлена в виде
Ьк+1 (Ук,...,Ук-г) = ^дз!к+1,з(Ук,...,Ук-г)+ 3+1
+ £
г=0
9Ък+\( Ук,..., Ук-г) ^ д д/к+1,з( Ук,..., Ук-г)
д Ук-г
3 = 1
дУ к-
к
или
Ък+\{ Ук,..., Ук-г) = Е 1=1 дз!к+1,з( Ук,..., Ук-г) + т+ъ
Я к+1 = Е 1=0
9Ьк+г(ук,...,Ук-г), 9ук-г
(18)
ЕП д 9(Ук,...,Ук-т) 3=1 Д 9ук-1
к-
(19)
где к = 0,1, 2,... — 1 и Ук-Г = 0 при к < г.
Для рассмотренных выше примеров разностные уравнения (18), описывающие результаты наблюдений, принимают частный вид [10].
1. Для модели, описывающей результаты наблюдений колебаний линейного осциллятора с затуханием
Ук = а0 вхр(—атк) еов(штк + фо) + ек,
в соответствии с разностными уравнениями (5) получаем
Уо = Дз + £о, У1 = + £1,
Ук = Д1Ук-1 + Д2 Ук-2 + Лк+1,
Щ+1 = — Д2 £к-2 — Д1 £к-1 + £к, к = 2, 3, 4,...,И — 1.
2. Для модели, описывающей результаты наблюдений колебаний системы с турбулентным трением
Ук = Л а° , со$(штк + фо) + £к, 1 + атк
в соответствии с разностными уравнениями (8) получаем
Уо = Дз + £о, У1 = + £1,
Ук + Ук-2 = Д1Ук-1 — Д2[кУк — Д1(к — 1) Ук-1+
+( к — 2)у к-2] + г)к+1,
щ+1 = [1 + (к — 2)Д2] £к-2—
—Д1[1 + (к — 1)Д2] £к-1 + (1 + к\2) £к, к = 2,...,Ы — 1.
3. Для модели, описывающей результаты наблюдений процессов с логистических трендом в форме функции Верхулста
(20)
Ук =
ао
+ £к,
1 + а1 ехр(—а1к) в соответствии с разностными уравнениями (16) получаем
' Уо = Дз + £о, —3 Уо + 4 У1 — У2 = Д1Уо + \2Уо + 42, Ук — Ук-2 = Д1 Ук-1 + Д2 у1-1 + Пк+1, т = —(3 + Д1 + 2Д2Уо) £о + 4£ 1 — £2, , Ш+1 = — £к-2 — ( Д1 + 2 Д2 Ук-1) £к-1 + £к, к = 2, 3,...,Ы — 1.
к = 2, 3,...,Ы — 1, (21)
На основе построенной системы разностных уравнений (18) и ее частных случаев (19)—(21) с помощью аппарата матричной алгебры формируется обобщенная регрессионная модель
b = FX + V = Р\£,
(22)
где Л = (А1, \2,..., Лга)т — вектор неизвестных коэффициентов системы разностных уравнений; е = (е0, £1, е2,..., £м-1)т — вектор разброса результатов наблюдений относительно модели (случайной помехи); ц = (щ, ц2, щ,..., Цм)т — вектор невязки (эквивалентного случайного возмущения);
ь = [bl(Уо), Ь2(Уо,Уl),. . ., br+1(Уо, У1, . . ., Уг), . . ., Ьм(Ум-г, Ум-r+1, . . ., Ум^ — вектор левой части обобщенной регрессионной модели; Р — регрессионная матрица размера [ N хп], элементы которой описываются формулами
!к,] = !к,](Ук-1, Ук-2,..., Ук-г-1), к = 1^, 2 = 1, п;
Р\ — матрица линейного преобразования вектора случайной помехи размера [ N хN], элементы которой описываются формулами
Pk,j
дЪк(Ук—1, ..., Ук-r-l)
9yj-i
£
г=1
дfk,i(Ук—1,..., Ук-r-l) д Уз-1 '
j = к - г, к (23)
для к ^ г + 1. При 1 ^ к ^ г в формуле (23) следует положить j = 1, к и считать, что переменные ук-г-1 = 0.
Представленные выше соотношения для вектора b и матриц F и Р\ для каждой исследуемой нелинейной зависимости ук = у(тк, а1, а2,...,ап) конкретизируются в процессе их формирования [11].
1. Для модели, описывающей результаты наблюдений колебаний линейного осциллятора с затуханием
Ук = а0 exp(-атк) сов(штк + ф0) + £к, в соответствии с разностными уравнениями (21) получаем
Р-
А =
. ^ Ум - -1)\ F
1 0 0
0 1 0
-\2 - 1
0- -\2 -\1
0 0
1
У2
0 0
Уо У1
10 01 00 00
У N-2 У N-3 0 0
0 0 0
0 0 0 0
0 0 - \2 -\1 1
0
0
2. Для модели, описывающей результаты наблюдений колебаний системы с турбулентным трением
к=
ао
■ со$(сотк + фо) + £к,
1 + атк
в соответствии с разностными уравнениями (22) получаем
^ =
Ь = (Уо,У\,Уо +У ьУ 1 + У2,..., Ум-з + Ум-О 00 00 У1 —(2У2 — Д1У1)
У2 —(3 Уз — Д12 У2 +У1)
т
10 01 00 00
Ум-2 — 1) Ум-1 — Х^И — 2) Ум-2 + ^ — 3) Ум-з] 0 0
Р:
Л =
1 0 0
0 1 0
1 — \1(1 + \2) 1 + 2X2
0 1 + Д2 — Д1(1 + 2Д2)
0 0 1 + 2 Д2
0 0 0
0 0 0
1 + 3 Д2 — Д1(1 + 3 Д2)
(25)
1 + (И — 3) Х2 — Х1[1 + (Ы — 2)Д2] 1 + (^ — 1)Д2
3. Для модели, описывающей результаты наблюдений процессов с логистических трендом в форме функции Верхулста
Ук
ао
+ к,
1 + а1 вхр(—а1к) в соответствии с разностными уравнениями (21) получаем
Ь= (Уо, —3Уо + 4У1 — У 2, У 2 — Уо, Уз — Ум-1 — Ум-з)
т
0 0 1
о о2 0
^ = 1 о12 0
2 У122 0
м-2 м2 - 2 0
Р:
А —
1
-(3 + \1 +2Х2У0) -1 0 0
0
4
-(Ах +2 Х2У1) -1 0
0 0 0 1
-(Ах + 2 Х2У3)
1
0 -1 1
-(Ах +2 А2У2) 1
(26)
0 0 0 0 0
-(Ах + 2 А2 ум-2) 1
5. Итерационная процедура уточнения среднеквадратичных оценок коэффициентов обобщенной регрессионной модели. Из первого уравнения обобщенной регрессионной модели (22) с учетом второго уравнения при невырожденной матрице Р\ можно получить
VA — СаА + е,
где = Р-Ь, Сх = Р-^ и е = Р-
Среднеквадратичные оценки коэффициентов обобщенной регрессионной модели, как отмечалось выше, находятся из условия минимизации остаточной суммы квадратов:
2
2
Г\\ — \\у-у\\ ^ mm.
С учетом преобразованного первого уравнения регрессии в (22) имеем
\к\\2 — \\^а - СаА\\2 — (г/А - СаА)т( и а - СаА) ^ min.
Отсюда получаем
ете — vJvA - 2иТСАА + АтС~ТСаА ^ min.
Используя аппарат матричной алгебры, можно показать справедливость следующих формул дифференцирования.
Пусть элементы матриц U(А) размера [Nхт] и V(А) размера [тхр] зависят от элементов вектора А размера [пх1]. Тогда имеет место формула
а[U w (x)]_au (Х)Вг м
дА
дА
дА
(27)
где
д [U (\)V (А)] ' dUV дUV\ \дUV'
дА дАх дА2\ \ дАп
блочная матрица-строка разме-
ра [Nхрп] производных от произведения матриц U^)V(А);
0
0
0
ди " ди ди \ \ ди '
дЛ = дЛ1 дЛ2\ \ дЛп
дУ ' дУ дУ ! \ дУ '
дЛ = дЛ1 дЛ2\ \ дЛп
(Л) =
у 0
0
0
у
0
у
— блочная матрица-строка размера [ Жхтл];
— блочная матрица-строка размера [ Nхрп];
— блочно-диагональная матрица размера [ тпхрп];
0 — нулевые матрицы размера [ тхр].
Рассмотрим производную по вектору Л € е € . Можно показать справедливость формулы
д( еТ е)
п от скалярной величины ете,
д Л
д
= 2£ т дЛ-
где результат дифференцирования есть матрица-строка размера [1хп]. С учетом равенства е = и\ — С\Л и формулы (29) отсюда получаем
д [(г/Л — СхЛ)т( — вхЛ)]
д Л
= 2( их — ОхЛ)
т д (^ — СлЛ)
д Л
2(- -™т(дЛ - ^Л
где Л =
Л 01 01
01 Л 01
01 01 Л
— блочно-диагональная матрица размера [п2хп],
0 — нулевые векторы размера [ N х1].
Так как система нормальных уравнений при среднеквадратичном оценивании в векторной форме обычно записывается в виде матрицы-столбца:
~д( ете) тт
д( Л) = qrad Q(Л) = получаем следующее равенство:
9( Л) = 2
д Л
(дVА д°Х Л т От
1\Ж — жЛ) — От
0,
(их — ОхЛ) = 0.
Отсюда система нормальных уравнений для задачи среднеквадратичного оценивания принимает вид
От (дд°\ Лт"
От \Ж — жЛ
°\Л =
От (дд°\ Лт
(28)
Рассмотрим два частных случая, вытекающих их этой формулы. Пусть и (Л) —квадратная невырожденная матрица размера [Ж хЖ ]. Тогда имеет место формула производной от обратной матрицы:
ди-1(Л)
= —и-1( Л)
ди (Л) ' дЛ
Ви-1 (Л),
(29) 687
где Bv-1 (X) =
'и-1 0 0
0 и-1 0
0 0 и-1
блочно-диагональная матрица разме-
В соответствии с принятым обозначением v\ = Р-1Ь и G\ = P-lF с уче-
ра [Nn х Nn], 0 — нулевые матрицы размера [N xN].
В соответствии с приняты] том формул (27) и (29) имеем
dvx dG\
dX dX
-Л
дРА" 1 дРА" 1 дРА" 1
Bb--^—Bp Л = ———B,q =
д X
д X
д X
= -P -1 ^BP-iB„ = -P-1 d-P±BF
д X
где B£(X) =
'Р-\Ъ -FX) 01 01
01 Р-\Ъ -FX) 01
01 01 Р-1 (Ь -FX)
д X
— блоч-
но-диагональная матрица размера [ Мпхп], 0 — нулевые векторы размера [Мх 1];
дРх \дРх\дРх\ \дРх 1 . г,, лг ,
—— = ———: •• ■ — —блочная матрица-строка размера хМп\. дХ дХ1 \ дХ2\ ■. дХп \
С учетом этого равенства имеем гт (д"\ дО\ .\т =
=рт(р-у+вт( дХ )т (р. у={р+^У^ У-
Тогда из (28) получаем систему нормальных уравнений для задачи среднеквадратичного оценивания коэффициентов обобщенной регрессионной модели (22):
(р + жВ^У^ 1рХ = (р + жВ^У^1 ь' (30)
где матрица Од находится по формуле Од = Р\Р\ .
Обобщением полученных формул (28) и (30), описывающих системы нормальных уравнений, является решение задачи среднеквадратичного оценивания коэффициентов обобщенной регрессионной модели общего вида
Ьх = F\X + г], V = Р\ £,
(31)
в которой элементы вектора левой части и матрицы регрессоров зависят от коэффициентов Х^ (например, как в матрице регрессоров (25) для систем с турбулентным трением). Для этого, полагая в (30) = Р^1 Ъ\ и С\ = = Р^ 1РЛ, с учетом формул (27) и (29) получаем
<9г/Л <9СЛ
дА дА
Л
дР-1 .л
+ дА -
-196 Л 9Р;
-1
дА
Л -Р:
-1
дА
Л
^+р-( ^ - дА л)
= —1К'* + £ - Ж Л)
-»-и др±п + дп
- Рх I дА * дА + дА
Л
" д Ъх дЬ х \ \д Ъх 1
дА1 дА2\ ' : дАп
где — - ^(А) -
для вектор-функции Ь\; строка размера [ Ж хп2]. Отсюда имеем
дРх дРх дРх\ \дРх
дА дА1 дА2\ : дАп
матрица Якоби размера [ N хп] — блочная матрица-
гт двх Лт ( +дрхя 9ЬХ + дрх т(р-1лт
-1 Ж- ЖЛ) -гл + ж* - ж + жЛ) {р-} •
Подставляя полученный результат в (28), получаем систему нормальных уравнений для задачи среднеквадратичного оценивания коэффициентов обобщенной регрессионной модели (31):
(* + ^ - Ж + дАЛ)Т^А -
- + Ж*'- - дА + ЖА)т°-1 ^ (32)
Формула (32) обобщает полученные выше результаты, описанные соотношениями (28) и (30). В частности, при Рх = Е — единичной матрице, - 0 — нулевой матрице, и формула (32) вырождается в соотношение (28). При 6 и Р, не зависящих от А^, ^ - 0 и ^^ - 0, а формула (32) обращается в формулу (30).
Если предположить, что элементы вектора Ь\, матриц Рх и Рх, вычисленные при некотором приближении Л(г) вектора коэффициентов, — константы и не зависят от А, то в (32) можно положить ^ - 0, ^^ - ^^ - 0. Тогда система нормальных уравнений для задачи среднеквадратичного оценивания коэффициентов обобщенной регрессионной модели принимает наиболее простой вид:
■%) РХИ)А - Р\а)ь\(о • (33)
Очевидно, что при заданных значениях вектора коэффициентов Л(г) система нормальных уравнений (33) линейна относительно вектора А, и ее решение при невырожденной матрице Р^а) ПГ^Р^фА, зависящее от А(г), может быть найдено по формуле
Л - (Рха) РХа)А) 1рТ) П-^ Ь-Ха) •
В общем случае система нормальных уравнений (32) является нелинейной. Ее решение может быть найдено методом простых итераций после преобразования к следующему виду:
Х =
1
(4+
дРх„ дЪх дРх чт
Жв*—Ж+ЖА) О- ^
Выбирая начальное приближение вектора коэффициентов обобщенной регрессионной модели, итерационную процедуру уточнения оценок Х(г) вектора коэффициентов можно записать в виде
Х(Ш) = (Щи П^^ХУнТ) О-г) ЬХ(г), 1 = 0,1,2,..., (34)
где матрица Н\, вычисленная в точке Х(г), описывается выражением
н\0) = +
д^Ш в® - ^ + ДО0, г = о, 1,2,... дХ дХ дХ
Очевидно, что для системы нормальных уравнений (33) выполняется Н—^ =
= ¿) и итерационная процедура уточнения оценок Л(г) вектора коэффициентов описывается формулой
А(т) = (рх(г) ^^^Г-) ^ Ъм, г = 0 12,...
А(^ г) О ,
В рассмотренных выше примерах построения обобщенной регрессионной модели колебаний линейного осциллятора с затуханием и модели, описывающей результаты наблюдений процессов с логистическим трендом в форме функции Верхулста, итерационная процедура уточнения оценок вектора коэффициентов принимает наиболее простой вид:
Х(+ = (р ТО^)РХ)-1РТО—) Ь, 1 = о, 1, 2,...
А( ^
(35)
где матрица Од(^ = Р^^Р-—) вычисляется по формулам (24) или (26) соответственно.
Проведенные численно-аналитические исследования показали, что аппроксимация нелинейной системы нормальных уравнений (32) линейной системой (33), элементы которой вычисляются в точке Л(г), с последующей итерационной процедурой (35) уточнения оценок вектора коэффициентов обеспечивает величину суммы квадратов отклонений модели от результатов наблюдений \\у — у\\2, которая, как правило, отличается от минимально возможного значения не более чем на 1-2 %. Для уточнения оценки Х(г+1) до величины не более 1 %:
\\х(<+1) — Х«\\ < о.о1 ■ \\\(*+1)\\
обычно требуется менее десяти итераций.
Оценку начального приближения Х(°) вектора коэффициентов системы разностных уравнений можно найти различными способами в зависимости от вида обобщенной регрессионной модели. Например, для обобщенной регрессионной модели вида (22) оценку начального приближения Х(°) достаточно
просто можно наити из условия минимизации невязки:
IHI2 = ||6-FA||2 ^ min
по формуле A(0) = (FTF)-lFTb.
При более сложной структуре разностных уравнении, например, когда элементы матрицы регрессоров F\ зависят от коэффициентов разностного уравнения, исходную регрессионную модель b = F\A + г] можно аппроксимировать моделью вида b = F*A* + rj за счет увеличения числа коэффициентов в модели. Например, для систем с турбулентным трением разностное уравнение
Ук + Ук-2 = Aiyk-i - A2 [кук - Ai(k - 1)ук-i + (k - 2)ук-2] + Vk+i,
которое содержит два несвязанных между собой коэффициента Ai и A2, можно представить в виде трехпараметрическои модели
Ук + Ук-2 = Aiyk-i - A2 [kyk + (k - 2)ук-2] + A3(k - 1)ук-i + щ+i. (36)
Очевидно, что разностное уравнение (36) содержит уже три коэффициента, причем на коэффициент A3 наложено ограничение вида A3 = AiA2. Аппроксимация исходного разностного уравнения моделью (36) без учета ограничений на коэффициент A3 позволяет использовать в качестве оценки начального приближения A(0 величину, найденную по формуле
A(0) = (F*TF *)-iF *т b,
где матрица F*, размер которой с учетом еще двух коэффициентов в (20) равен [ N х 5], формируется в соответствии с формулой (36).
6. Вычисление оценок параметров нелинейной математической модели на основе среднеквадратичных оценок коэффициентов разностных уравнений. В основе вычисления оценок параметров нелинейной зависимости у = f(t, ai,a2,..., ап) лежат соотношения (3), полученные при построении системы разностных уравнений (2). Изначально разностные уравнения строятся таким образом, чтобы нелинейная система (3) имела простой вид и легко разрешалась относительно искомых параметров.
В рассматриваемом выше примере нелинейной математической модели (4) линейного осциллятора с затуханием система нелинейных уравнений принимает вид (6). Из нее легко можно получить следующие соотношения для вычисления оценок параметров нелинейной зависимости (4):
- i w ^ - i Ai
а =--ln(—A2), ш = -arceos— „ ,
2 2 - A2
A
/\2 ( AiA3 - 2A4)2 Г i AiA3 - 2 A4
ao = \ ¡A2--^——^—, фо = arctg-
4 A2 + Ai -A -Ai
Для модели, описывающей ординаты свободных колебаний диссипатив-ной системы с турбулентным трением, с учетом формул (9) получаем соотношения для вычисления оценок параметров нелинейно зависимости (7):
. Х" „ 1 Х1 л
а = —, ш = — агсссе —, а° = л т т 2 \
Х [Х1Х3 — 2 Х4(1 + Л2)]:
Х3 +---
3 4 — Х1
Х Х1Х3 — 2Х4(1 + Х2)
Хзу/ 4 — Х1
И, наконец, для модели, описывающей результаты наблюдений процессов с логистическим трендом в форме функции Верхулста, в соответствии с формулами (15) соотношения для вычисления оценок параметров нелинейной зависимости (13) принимают вид
Х1 Х1 Х1 ° Х1
• а° = — Х2 • а1 =—еП -2°) + ХХТз) ■
7. Оценка погрешности результатов вычислений на основе методов статистической обработки данных эксперимента. Оценка погрешности результатов вычислений параметров математической модели является важнейшим и обязательным шагом в алгоритме численного метода на основе разностных уравнений. Без этого заключительного этапа задача построения математической модели не может считаться решенной.
Так как разностные уравнения, описывающие результаты наблюдений, содержат ненаблюдаемые переменные , описывающие естественный разброс результатов эксперимента у к относительно построенной модели у к, величины ек = Ук — Ук можно рассматривать как случайную аддитивную помеху (возмущение). При таком подходе оценку погрешности результатов вычислений следует проводить на основе известных статистических методов обработки экспериментальных данных [1-3]. В первом приближении будем предполагать, что математическая модель в форме построенных разностных уравнений и соответствующая им обобщенная регрессионная модель априори адекватны наблюдаемому процессу, причем случайные величины к имеют нормальный закон распределения, нулевые математические ожидания, одинаковые дисперсии а"2 и некоррелированны между собой, то есть матрица дисперсий-ковариаций вектора случайного возмущения имеет вид
V [е] = Еа\
где Е — единичная матрица [1]. С учетом принятых допущений предлагается следующий алгоритм процедуры оценивания погрешности результатов вычислений параметров математической модели в форме нелинейной функциональной зависимости:
- вычисление оценки дисперсии 2[ ] случайной помехи в результатах наблюдений;
- построение матрицы дисперсий-ковариаций V[Х] вектора оценок коэффициентов обобщенной регрессионной модели;
- построение матрицы дисперсий-ковариаций V[а] вектора оценок параметров нелинейной зависимости;
- вычисление доверительных интервалов для оценок параметров нелинейной зависимости;
- вычисление оценок дисперсии з2[Х] результатов вычислений на основе построенной математической модели;
- построение доверительных интервалов для оценок результатов вычислений у к на основе построенной математической модели.
Оценка в2 [ек] дисперсии а2 случайной помехи в результатах наблюдений может быть найдена различными способами в зависимости от априорной информации. В частности, на основе остаточной суммы квадратов Qres = ете, где е = Р-1 Ьд — Р-1Р^Л — остатки, служащие оценками случайной величины
е, за оценку дисперсии в2[£к] может быть принята остаточная дисперсия [1]: 8 2[ек] = 8 1е8, где
2 = =( ъ% — рхХ)тп-1( Ъ% — рх\)
в гез N — п N — п '
а Х — вектор оценок коэффициентов, вычисленных на основе предлагаемого численного метода; V[Х] —объем выборки результатов наблюдений; N — число коэффициентов в обобщенной регрессионной модели.
При построении матрицы дисперсий-ковариаций V[Л] вектора оценок коэффициентов обобщенной регрессионной модели воспользуемся полученным соотношением (34), которое при завершении итерационной процедуры уточнения оценок вектора коэффициентов г ^ те принимает вид
- , т 1 чт т 1 дР\ -Ь\ 9Р\ , ч
Х = ^Г^ХЧ-1 Я = р-х + — + А.. (37)
Так как Х * Л, из (37) можно получить
Л * (ЯЛтП-1Рл)-1ЯЛтП-1 ЪХ = (Ялт0-1Рл)-1ЯлтП-1(РлЛ + г?) =
= (ЯлтП-1Рл)-1Ялт П--1РХЛ + (ЯлтО-1РА)-1Ялт(Рл-1)тРл-1^ =
= Л + (ЯлХ1Рл)-1Ялт(Рл-1)те.
Поскольку матрица (Я'^Р\)-1 Я'^(Р-1 )т не зависит от случайной величины е, при принятом выше допущении М[е] = 0, где М[ ■ ] —оператор математического ожидания:
М [Л] = Л + (ЯЛтО-1Рл)-1ЯЛт(Р-1)тМ И = Л,
то есть оценка ЛХ — несмещенная оценка вектора коэффициентов.
Матрицу дисперсий-ковариаций V[Л] находим по следующей формуле [1]:
V[Л] = М [(Х — М[Х])(Л — МЙ)т] = М[(Л — Л)(Х — Л)т] =
= М [(ЯЛтО-1РЛ)-1ЯЛт(РЛ-1)т^тРЛ-1ЯЛ(ЯЛтО-1] = = (ЯлтО-1Рл)-1Ялт(Р-1)тМ [^т]Рл-1Ял(ЯлтОл)-1 = = (ЯлтО-1Рл)-1Ялт(Р-1 )тРл-1Ял(Ялт Ол)-1а2 = (ЯлтО-1Рл)-1а£2.
Отсюда получаем формулу для вычисления матрицы оценок дисперсий-ковариаций результатов вычисления коэффициентов обобщенной ререссион-ной модели:
F [Л] = (HjQ-'Fx)-1 s2res.
Диагональные элементы матрицы V[Л] есть оценки дисперсий коэффици-
ентов \j :
v33 = s Г[Лз Ь
а недиагональные — оценки ковариаций между Хг и :
= ссу[Хг, Хз].
Пусть ак(Х1, Х2,..., Хп) — один из параметров нелинейной модели, известным образом зависящий от коэффициентов системы разностных уравнений Х = (Х1, Х2,..., Хп)т, а Хк(Хь Х2,..., Хп) — несмещенная оценка этого параметра ( М[Хк] = ак). Очевидно, что отклонение оценки от ее математического ожидания с учетом несмещенности оценок коэффициентов М [Хг] = Х1 можно представить в виде
п да п да
Хк — м[ак] = Хк — ак дХ(^ — Х*) = ^2 дХк(^ — мХ])-
1=1 г 1=1 г
Рассмотрим ковариацию между оценками г^к и а с учетом их несмещенности:
ссу[ак, а ] = м [(ак — м [Хк })(а1 — м [ак = м [(Хк — ак )(а1 — аг)] =
м [ Е dak& - мЕ длЛ - мЛ]\
г=1 ' 'г 3=1 ' 3
dak л -м^ ^dai
dak dai дЛг дЛз
м[ЕЕдЛгщЛ -м- мл])
dakdai1,/t и Л „win xr
d Л з
=1 =1
=1
Отсюда
1 п я я
+ 2 ЕЕ дЛ Ж-м [Л -м [А,])(Аз -м [Лз ])
г=1 з=г+1 г 3
' 1 '
сс[акМ = £ ^ ^а2[.Л+2 £ £ ^ ^ (38)
г=1 J г=13=1+1
В частном случае, когда I = к, имеем формулу для вычисления дисперсии оценки параметра Хк:
![ак ] = ^
г=1
/ дак V
™-1 п я я ![А.]+2 £ £ ^ ^ соу[А(, А,- ]. (39)
г=1 ]=+1
Рассмотрим теперь матрицу Якоби производных элементов вектора а = = ( а1, а2,..., ап)т по коэффициентам обобщенной регрессионной модели — элементам вектора А = (А1, Л2,..., Ап)т:
^а(А) =
д а1 д а1 д а1
дА1 да2 дА2 да2 дАп да2
дА1 дА2 дАп
дап дап дап
.дА1 дА2 дАп
Пусть "к и "г — два столбца этой матрицы:
"к =
/дак дак 1дА1' дА2'
дак\т дАпУ '
" =
/ даг даг
1дА7 'дА2'
даг \т дАпУ '
Тогда формулы (38) и (39) в матричной форме принимают следующий вид:
еоу[ак, аг ] = ад^ТУ [Аод ], &2[ак ] = "У [А" ].
Объединяя эти два соотношения, получаем формулу, описывающую матрицу оценок дисперсий-ковариаций результатов вычисления параметров нелинейной функциональной зависимости:
V [а] = WаV [А]^аТ
При построении доверительных границ погрешности результатов вычислений оценок параметров можно воспользоваться формулой
Дак = Ь/зв [ак ],
где «[ак] —оценка среднеквадратичного отклонения результата вычислений параметра ак. В первом приближении можно считать, что статистика Ьр имеет распределение Стьюдента с V = N — п степенями свободы. В этом случае при доверительной вероятности р = 0.95 и числах степеней свободы V ^ 25 можно принять Ьр = 2.1. Так как при нелинейной зависимости функция распределения погрешности результата вычислений неизвестна, доверительный интервал для случайной погрешности целесообразно строить не на основе распределения Стьюдента, а с использованием неравенства Чебышева [7]. Считая распределение погрешности результатов вычисления симметричным и одномодальным, неравенство Чебышева при построении доверительного интервала для оценки параметра а к запишем в следующем виде [7]:
4
Р(|Дак| ^ ¿ра[ак]) < -2 .
9
3
В этом случае границы доверительного интервала вычисляются при величине £3, которая находится из формулы доверительной вероятности р = 1 — 92.
В частности, при ß = 0.95 величина tß = 2.98. Найденная по формуле Aäk = = tßS [äk] величина может рассматриваться в качестве предельной абсолютной погрешности (с заданной вероятностью) результата вычислений параметра ak:
ak = äk ± Aaak.
Оценку предельной относительной погрешности можно найти по формуле
öäk = Aäk /\a k \.
При вычислении оценок дисперсии s2[yk] результатов вычислений на основе построенной математической модели у = y(t, ä\,ä2,..., än) можно воспользоваться соотношением, аналогичным приведенному выше:
s 2[yk ] = wJV [ä]wk,
где Wk — вектор-столбец частных производных нелинейной функциональной зависимости по параметрам, вычисленным в точке k:
Wk = G^,---,^)" * = 0,1, 2,...,N - 1. \oai oa2 oanJ
Если ввести матрицу Якоби Wy размера [N х п], элементы которой описываются формулой
wkj = т-1, * = 0,1, 2,... ,N - 1, j = 1, 2,... ,п,
Odj
то матрицу оценок дисперсий-ковариаций результатов вычислений ä на основе построенной модели у = y(t, ä\,ä2,..., än) можно представить в виде
V [ä] = WyV [aW. (40)
Диагональные элементы иц = w"V[a]wk матрицы (40) есть оценки дисперсии s2[yk] результатов вычислений yk на основе построенной модели.
Доверительный интервал для результата вычисления yk находится по формуле
A k = ß [ k],
где величина ß выбирается в зависимости от доверительной вероятности ß либо на основе распределения Стьюдента, либо на основе неравенства Чебы-шева, как описано выше.
Заключение. В данной работе описаны теоретические основы и алгоритм нового численного метода нелинейного оценивания, в основе которого лежат разностные уравнения, описывающие результаты наблюдений. Известные соотношения, связывающие коэффициенты обобщенной регрессионной модели с параметрами нелинейной функциональной зависимости, позволяют свести задачу среднеквадратичного оценивания параметров нелинейной модели к итерационной процедуре уточнения решения линейной регрессионной задачи. Предлагаются различные подходы к построению систем разностных уравнений, описывающих результаты наблюдений; построены новые соотношения, лежащие в основе итерационной процедуры уточнения оценок коэффициентов обобщенной регрессионной модели и позволяющие обеспечить минимум суммы квадратов отклонения результатов наблюдений от результатов
расчета по модели. Описанная в работе методика оценивания погрешности результатов вычислений обеспечивает завершенность и полноту процедуры параметрической идентификации, что позволяет рекомендовать разработанный численный метод для решения практических задач математического моделирования в различных областях научной и производственной деятельности.
Разработанный численный метод нелинейного оценивания на основе разностных уравнений имеет ряд преимуществ по сравнению с известными методами нелинейной регрессии [2,3]. Во-первых, при его применении не требуется нахождения производных первого и второго порядков от нелинейной функции, описывающей математическую модель. Во-вторых, одна из важнейших проблем нелинейного оценивания — выбор начального приближения, обеспечивающего сходимость итерационной процедуры уточнения оценок параметров модели — в предлагаемом методе успешно решается на основе минимизации невязки в первом линейном уравнении обобщенной регрессионной модели.
Для применения разработанного численного метода на основе разностных уравнений, как правило, достаточно лишь сформировать вектор , матрицу регрессоров F и матрицу Р\ линейного преобразования вектора случайной помехи в обобщенной регрессионной модели (22). Этот вектор и матрицы построены и приведены в целом ряде работ, посвященных применению описанного численного метода в задачах параметрической идентификации математических моделей процессов различной физической природы, в частности: линейных динамических систем [11,18,19], нелинейного дифференциального оператора второго порядка [16], математических моделей в форме нелинейных дробно-рациональных зависимостей, нелинейных диссипативных механических систем [11], в том числе систем с турбулентным трением [13]. Результаты апробации метода нелинейного оценивания на основе разностных уравнений также можно найти в работах [20-22], посвященных проблеме построения по результатам эксперимента моделей ползучести разупрочненно-го материала в пределах всех трех стадий ползучести. В этих же работах можно найти результаты численно-аналитических исследований и сравнительный анализ с известными методами нелинейного оценивания. Результаты апробации и численно-аналитических исследований на основе компьютерного моделирования подтверждают новизну, достоверность и практическую значимость представленных в данной работе математических моделей в форме разностных уравнений и описанных алгоритмов численного метода.
Конкурирующие интересы. Конкурирующих интересов не имею. Авторский вклад и ответственность. Я несу полную ответственность за предоставление окончательной версии рукописи в печать. Окончательная версия рукописи мною одобрена.
Финансирование. Работа выполнена при поддержке Российского фонда фундаментальных исследований (проект № 16-01-00249_a).
Библиографический список
1. Вучков И., Бояджиева Л., Солаков О. Прикладной линейный регрессионный анализ.
М.: Финансы и статистика, 1987. 238 с.
2. Draper N. R., Smith H. Applied regression analysis/ Wiley Series in Probability and Mathematical Statistics. New York etc.: John Wiley & Sons, 1981. xiv+709 pp. doi:
9781118625590.
3. Демиденко Е. З. Линейная и нелинейная регрессии. М.: Финансы и статистика, 1981. 302 с.
4. Bjorck A. Numerical methods in matrix computations / Texts in Applied Mathematics. vol.59. Cham: Springer, 2015. xvi+800 pp. doi: 10.1007/978-3-319-05089-8.
5. Bard Y. Nonlinear parameter estimation. New York: Academic Press, 1974. x+341 pp.
6. Gunst R. F., Mason R. L. Regression analysis and its application. A data-oriented approach / Statistics: Textbooks and Monographs. vol. 34. New York, Basel: Marcel Dekker, 1980. xv+402 pp.
7. Грановский В. А., Сирая Т. Н. Методы обработки экспериментальных данных при измерениях. Л.: Энергоатомиздат, 1990. 288 с.
8. Marquardt D. W. An algorithm for least-squares estimation of nonlinear parameters // J. Soc. Ind. Appl. Math., 1963. vol.11, no. 2. pp. 431-441. doi: 10.1137/0111030.
9. Hartley H. O., Booker A. Nonlinear least squares estimation// Ann. Math. Stat., 1965. vol.36. pp. 638-650. doi: 10.1214/aoms/1177700171.
10. Формалиев В. Ф., Ревизников Д. Л. Численные методы. М.: Физматлит, 2006. 400 с.
11. Зотеев В. Е. Параметрическая идентификация диссипативных механических систем на основе 'разностных уравнений. М.: Машиностроение, 2009. 344 с.
12. Деч Г. Руководство к практическому применению преобразования Лапласа и z-преобразования. М.: Наука, 1971. 288 с.
13. Егорова А. А Метод параметрической идентификации систем с турбулентным трением / Математическое моделирование и краевые задачи: Труды Восьмой Всероссийской научной конференции с международным участием. Часть 2. Самара: СамГТУ, 2011. С. 143-156.
14. Волков Е. А. Численные методы. СПб.: Лань, 2004. 256 с.
15. Дубовцев А. В., Ермолаев М. Б. Прогнозирование развития рынка мобильной связи на основе 5-образных моделей // Современные наукоемкие технологии. Региональное приложение, 2010. №4 (24). С. 39-41.
16. Martino J. P. Technological forecasting for decisionmaking. New York: American Elsevier, 1972. xviii+750 pp.
17. Дуброва Т. А. Статистические методы прогнозирования. М.: Юнити-Дана, 2003. 206 с.
18. Зотеев В. Е. Параметрическая идентификация линейной динамической системы на основе стохастических разностных уравнений // Матем. моделирование, 2008. Т. 20, №9. С. 120-128.
19. Зотеев В. Е., Стукалова Е. Д., Башкинова Е. В. Численный метод оценки параметров нелинейного дифференциального оператора второго порядка // Вестн. Сам. гос. техн. ун-та. Сер. Физ.-мат. науки, 2017. Т. 21, №3. С. 556-580. doi: 10.14498/vsgtu1560.
20. Зотеев В. Е., Макаров Р. Ю. Численный метод определения параметров модели ползучести в пределах первых двух стадий // Вестник Самарского университета. Аэрокосмическая техника, технологии и машиностроение, 2017. Т. 16, №2. С. 145-146. doi:10.18287/2541-7533-2017-16-2-145-156.
21. Макаров Р. Ю. Численный метод определения параметров кривой ползучести на основе закона Содерберга // Вестник Самарского университета. Аэрокосмическая техника, технологии и машиностроение, 2015. Т. 14, №2. С. 113-117. doi: 10.18287/ 2412-7329-2015-14-2-113-118.
22. Зотеев В. Е., Макаров Р. Ю. Численный метод оценки параметров деформации ползучести при степенной зависимости параметра разупрочнения // Современные технологии. Системный анализ. Моделирование, 2016. Т. 51, №3. С. 18-25.
Vestn. Samar. Gos. Tekhn. Univ., Ser. Fiz.-Mat. Nauki
[J. Samara State Tech. Univ., Ser. Phys. Math. Sci.], 2018, vol. 22, no. 4, pp. 669-701
d http://doi.org/10.14498/vsgtu1643
ISSN: 2310-7081 (online), 1991-8615 (print)
MSC: 65C20, 65P40
A numerical method of nonlinear estimation based on difference equations
V. E. Zoteev
Samara State Technical University,
244, Molodogvardeyskaya st., Samara, 443100, Russian Federation.
Abstract
The article considers a new numerical method for estimating the parameters of nonlinear mathematical models based on difference equations describing the results of observations. The algorithm of the numerical method includes:
- the construction of a linear-parametric discrete model of the process under study in the form of difference equations, the coefficients of which are known to be associated with the parameters of a nonlinear mathematical model;
- the formation of a generalized regression model based on the difference equations;
- the calculation of the initial approximation estimate and the iterative procedure for refining the mean-square estimates of the coefficients of the generalized regression model;
- the calculation of the estimates of the parameters of the nonlinear mathematical model based on the mean-square estimates of the coefficients of the difference equations;
- evaluation of the error of the results of calculations based on the methods of statistical processing of experimental data.
Various approaches to the construction of systems of difference equations for mathematical models in the form of nonlinear functional dependencies are proposed. The relations underlying the iterative process of refining the coefficients of the generalized regression model constructed on the basis of difference equations are obtained. The procedure for estimating the error of the results of calculations of the parameters of nonlinear functional dependencies, which are known to be associated with the coefficients of the system of difference equations, is described. The application of the numerical method based on the difference equations is illustrated by the examples of estimation of the parameters of the mathematical model of the linear oscillator with attenuation, the model of free oscillations of the dissipative mechanical system with turbulent friction, as well as the parameters of the logistic trend described by the Verhulst (Pearl-Reed) function.
Research Article
3 ©® The content is published under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/) Please cite this article in press as:
Zoteev V. E. A numerical method of nonlinear estimation based on difference equations, Vestn. Samar. Gos. Tekhn. Univ., Ser. Fiz.-Mat. Nauki [J. Samara State Tech. Univ., Ser. Phys. Math. Sci.], 2018, vol. 22, no. 4, pp. 669-701. doi: 10.14498/vsgtu1643 (In Russian). Author's Details:
Vladimir E. Zoteev http://orcid.org/0000-0001-7114-4894
Dr. Tech. Sci.; Professor; Dept. of Applied Mathematics and Computer Science;
e-mail: [email protected]
Keywords: mathematical model, nonlinear regression analysis, system of difference equations, generalized regression model, mean square estimation, statistical processing of experimental results.
Received: 14th September, 2018 / Revised: 5th November, 2018 / Accepted: 12th November, 2018 / First online: 18th December, 2018
Competing interests. I declare that I have no competing interests.
Author's Responsibilities. I take full responsibility for submitting the final manuscript
in print. I approved the final version of the manuscript.
Funding. This work was supported by the Russian Foundation for Basic Research (project no. 16-01-00249_a).
References
1. Vuchkov I., Boyadzhieva L., Solakov O. Prikladnoi lineinyi regressionnyi analiz [Applied linear regression analysis]. Moscow, Finansy i statistika, 1987, 238 pp. (In Russian)
2. Draper N. R., Smith H. Applied regression analysis, Wiley Series in Probability and Mathematical Statistics. New York etc., John Wiley & Sons, 1981, xiv+709 pp. doi: 9781118625590.
3. Demidenko E. Z. Lineinaia i nelineinaia regressii [Linear and Nonlinear Regression]. Moscow, Finansy i statistika, 1981, 302 pp. (In Russian)
4. Bjorck A. Numerical methods in matrix computations, Texts in Applied Mathematics, vol.59. Cham, Springer, 2015, xvi+800 pp. doi: 10.1007/978-3-319-05089-8.
5. Bard Y. Nonlinear parameter estimation. New York, Academic Press, 1974, x+341 pp.
6. Gunst R. F., Mason R. L. Regression analysis and its application. A data-oriented approach, Statistics: Textbooks and Monographs, vol. 34. New York, Basel, Marcel Dekker, 1980, xv+402 pp.
7. Granovskii V. A., Siraya T. N. Metody obrabotki eksperimental'nykh dannykh pri izmereni-iakh [Methods of Processing of the Measurements of Experimental Data]. Leningrad, Ener-goatomizdat, 1990, 288 pp. (In Russian)
8. Marquardt D. W. An algorithm for least-squares estimation of nonlinear parameters, J. Soc. Ind. Appl. Math., 1963, vol.11, no. 2, pp. 431-441. doi: 10.1137/0111030.
9. Hartley H. O., Booker A. Nonlinear least squares estimation, Ann. Math. Stat., 1965, vol. 36, pp. 638-650. doi: 10.1214/aoms/1177700171.
10. Formaliev V. F., Reviznikov D. L. Chislennye metody [Numerical methods]. Moscow, Fiz-matlit, 2006, 400 pp. (In Russian)
11. Zoteev V. E. Parametricheskaia identifikatsiia dissipativnykh mekhanicheskikh system na osnove raznostnykh uravnenii [Parametric identification of dissipative mechanical systems based on difference equations]. Moscow, Mashinostroenie, 2009, 344 pp. (In Russian)
12. Dech G. Rukovodstvo k prakticheskomu primeneniiu preobrazovaniia Laplasa i z-preobrazovaniia [Guide to the practical application of Laplace transform and ^-transform]. Moscow, Nauka, 1971, 288 pp. (In Russian)
13. Egorova A. A. Method of parametric identification of the systems with turbulent friction, In: Mathematical Modeling and Boundary-Value Problems, Proceeding of the Eigth all-Russian Scientific Conference with International Participation. Part 2. Samara, Samara State Technical Univ., 2011, pp. 143-156 (In Russian).
14. Volkov E. A. Chislennye metody [Numerical Methods]. St. Petersburg, Lan', 2004, 256 pp. (In Russian)
15. Dubovtsev A. V., Ermolaev M. B. The forecasting of mobile communication market development on the basis of 5-shaped models, Sovremennye naukoemkie tekhnologii. Regionalnoe prilogenie, 2010, no. 4(24), pp. 39-41 (In Russian).
16. Martino J. P. Technological forecasting for decisionmaking. New York, American Elsevier, 1972, xviii+750 pp.
17. Dubrova T. A. Statisticheskie metody prognozirovaniia. Moscow, Yuniti-Dana, 2003, 206 pp. (In Russian)
18. Zoteev V. E. Parametrical identification of linear dynamical system on the basis of stochastic difference equations, Matem. Mod., 2008, vol.20, no. 9, pp. 120-128 (In Russian).
19. Zoteev V. E., Stukalova E. D., Bashkinova E. V. Numerical method of estimation of parameters of the nonlinear differential operator of the second order, Vestn. Samar. Gos. Tekhn. Univ., Ser. Fiz.-Mat. Nauki [J. Samara State Tech. Univ., Ser. Phys. Math. Sci.], 2017, vol.21, no. 3, pp. 556-580 (In Russian). doi: 10.14498/vsgtu1560.
20. Zoteev V. E., Makarov R. Yu. Numerical method of determining creep model parameters within the first two stages of creep, Vestnik of Samara University. Aerospace and Mechanical Engineering, 2017, vol.16, no. 2, pp. 145-146 (In Russian). doi: 10.18287/ 2541-7533-2017-16-2-145-156.
21. Makarov R. Yu. Numerical method for determining the parameters of a creep curve on the basis of Soderberg law, Vestnik of the Samara State Aerospace University, 2015, vol. 14, no. 2, pp. 113-117 (In Russian). doi: 10.18287/2412-7329-2015-14-2-113-118.
22. Zoteev V. E., Makarov R. Yu. Numerical method of estimation of parameters of deformation of creep in the exponential dependency of parametr of weakening from the strain, Modern technologies. System analysis. Modeling, 2016, vol.51, no. 3, pp. 18-25 (In Russian).
