Исследование эффективности применения линейных дискретных моделей при определении параметров математических моделей в форме обыкновенных дифференциальных уравнений Текст научной статьи по специальности «Математика»

УДК 534.282

В.Е.Зотеев

ИССЛЕДОВАНИЕ ЭФФЕКТИВНОСТИ ПРИМЕНЕНИЯ ЛИНЕЙНЫХ ДИСКРЕТНЫХ МОДЕЛЕЙ ПРИ ОПРЕДЕЛЕНИИ ПАРАМЕТРОВ МАТЕМАТИЧЕСКИХ МОДЕЛЕЙ В ФОРМЕ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ

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

Постановка задачи. При изучении явлений и процессов в природе, решении прикладных задач в самых разных областях науки и техники широко используются математические модели изучаемого объекта [1]. Они позволяют аналитически исследовать динамическую систему, определять ее оптимальные параметры, прогнозировать развитие процесса. При этом математическое моделирование предполагает использование современных средств вычислений с их развитым программным обеспечением.

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

Задача определения коэффициентов ОДУ (параметров объекта исследования) является обратной к задаче решения дифференциального уравнения. На практике она решается на основе обработки выборки результатов измерений сигнала на выходе динамической системы.

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

Физические процессы и их математическое описание. Рассмотрим некоторые процессы различной природы, описывающиеся одной и той же математической моделью в форме ОДУ [1].

Одна из простейших математических моделей эволюции биологической популяции осно-

дх

вана на предположении, что численность х популяции и скорость — изменения численности

дХ дх

во времени Х пропорциональны, т.е. — = кх . Коэффициент к можно трактовать как разность

Ж

коэффициентов рождаемости и смертности. Если считать коэффициент к постоянным, то решением этого уравнения будет зависимость х(х) = х0екХ, где х0 - численность популяции в момент времени Х = 0 .

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

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

ленности. В этом случае дифференциальное уравнение будет иметь вид

йх , — = к йі

ґ \ 1 --Х.

V

х„

х.

(1)

/

Рассмотрим математическое описание некоторых химических процессов [1]. Скорость химической реакции характеризуется скоростью образования конечного продукта реакции. Согласно закону действующих масс эта скорость при постоянной температуре пропорциональна произведению текущих концентраций реагирующих веществ. Пусть для получения х молекул

х

конечного продукта X требуется х молекул вещества А и молекул вещества В, т.е. уравнение

реакции можно записать в виде 2А+В=2Х. В момент времени Х = 0 начальные концентрации веществ А и В равны соответственно а и Ь, а концентрация конечного продукта х(0) = 0. В силу закона действующих масс в этот момент времени начальная скорость реакции пропорциональна произведению аЬ, а в некоторый момент времени Х > 0 при концентрации х(х) конечного продукта получим ОДУ

(2)

где к > 0 - константа скорости рассматриваемой реакции в данных условиях.

Математические модели (1) и (2) описывают совершенно различные по природе динамические процессы, однако имеют одинаковую структуру. В дальнейшем будем рассматривать ОДУ первого порядка вида

йх

йі

= а0 + а1 х + а2 х

(3)

■ 4а0а2 > 0.

которое обобщает ОДУ (1) и (2) с учетом условия В = а12

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

Іі ) = :

■ х2Се

4Ьі

1 - Се

(4)

где х1 и х2 - корни выражения, стоящего в правой части уравнения (3); С = -

произ-

вольная постоянная, зависящая от начальных условий х(0) = х0 .

Рассмотрим динамические процессы, которые описываются ОДУ второго порядка. В качестве примера остановимся на системе с гидравлическим демпфером, который создает сопротивление движению поршня, зависящее от скорости и пропорциональное ее первой степени (рис. 1).

_____-Ьх'(і) //////////////////////////////

< - сх(і)

•(і)

777777777777777777777777777777

Р и с. 1. Функциональная схема гидравлического демпфера

х

т

х0 х1

х0 - х2

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

Баланс сил в данной системе приводит к дифференциальному уравнению

тх"(ґ ) = —сх(ґ)- Ьх \ґ), (5)

где — сх(ґ) - восстанавливающая сила; — Ьх' (ґ) - диссипативная сила (вязкое трение). При малых значениях коэффициента Ь решение этого ОДУ представляет собой затухающие по экспо-

[С

ненте гармонические колебания. Обозначив через СО 0 =л— собственную частоту колебаний,

V т

а через 5 = — показатель затухания колебаний, получаем

2т

х" + 2 5% + Сс02 x = 0. (6)

Решение уравнения (6) может быть представлено в виде

х(ґ) = а0е~5ґ со$(а>ґ + у0), (7)

где (О = д/СО2 — 52 » (О0 - частота колебаний системы; а0 и у0 - амплитуда и фаза колебаний,

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

рЬ

ОДУ (5) по формуле 8 = . Параметры о и 8 являются динамическими характеристика-

Ыст

ми (ДХ) механической системы и часто используются в качестве информативных признаков ее технического состояния.

Разностные уравнения динамических процессов и алгоритм вычисления коэффициентов ОДУ на их основе. Разностные уравнения, связывающие результаты измерений наблюдаемого процесса с шагом (периодом дискретизации) %, строятся на основе ОДУ [2]. Производные заменяются по формулам численного дифференцирования х ( ) = ^к--------------—- , которые

2%

имеют второй порядок точности относительно % . При такой замене в ОДУ (3) получаем разностное уравнение вида

Хк *к—2 = а0 + а1 хк—1 + а2х1—1, к = 2,3,4... (8)

2%

Оно легко преобразуется в линейное относительно коэффициентов Яу , у = 1,3 уравнение

1 хк—1 + Я2хк—1 + Я3 = хк — хк—2 , к = 2,3,4... . (9)

Коэффициенты Я у связаны с коэффициентами а0, а 1 и а2 соотношениями

Я Я я2 (10)

а0 = —, а1 = —, а2 = —. (10)

0 2% 1 2% 2 2%

При замене производных в ОДУ второго порядка (5) получаем разностное уравнение вида

х — 2 хк—1 + х—2 + AХL-Хk-l + £хк—1 = 0, к = 2,3,4... .

% т 2% т

Отсюда после простых преобразований получаем

1 хк—1 + Я2хк—2 = хк, к = 2,3А.. . (11)

Коэффициенты разностного уравнения (11) связаны с коэффициентами ОДУ соотношениями

Ь = 2(1 + Я2) с = 2(1 — Я1 — Я2)

т % (1 — Я2) ’ т % (1 — Я2)

(12)

При этом динамические характеристики системы определяются по формулам

Л _ 2р(1 + Я2) Л1 д/2(1 -1 -Я>)-(1 + Я,)2

рт (1 - Я - Я2 )(1 - Я2) ’ т( - Я2) '

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

В данной работе решается обратная задача: по известным значениям функции х(х) при Xк _ кт , к _ 0,1,2,...N -1, где N - объем выборки, вычисляются коэффициенты ОДУ. Полагая в (9) или (11) к _ 2,3, ..^, получаем систему из N-1 линейных алгебраических уравнений (СЛАУ), содержащую в первом случае три, а во втором - две неизвестных.

Отсчеты наблюдаемого динамического процесса хк _ х(тк) всегда содержат некоторую погрешность, связанную, например, с точностью средств измерений. Поэтому при решении полученной СЛАУ используют широко распространенные в практике обработки экспериментальных данных методы прикладного линейного регрессионного анализа, например метод наименьших квадратов (МНК) [3].

Минимизация по переменным Яу сформированных в соответствии с МНК функционалов

для разностного уравнения (9) приводит к нормальной системе уравнений

N N N N

Я1 ^ ХАт-1 + Я2 ^ Хк-1 + Я3 ^ Хк-1 _ ^ (Хк - Хк-2 )Хк-1

к_2 к_2 к_2 к_2

N N NN

Яі Е Xk-i + Я2 Е Xk-i + Я3 Е Xk-i =Е (xk Xk-2 К-i , (14)

k=2 k=2 k=2 k=2

NN N

Яі Е Xk-і + Я2 Е Xk-i + Я3 ( - і) = Е (

k=2 k=2 k=2

а для разностного уравнения (11) - к нормальной системе уравнений

N N N

Яі Е Xk-i + Я2 Е Xk-2Xk-■ = Е XkXk-і >

k=2 k=2 k = 2

N NN

Яі Е Xk-■ Xk-2 + Я2 Е Xk-2 =Е XkXk-2 • (ІЗ)

k = 2 k=2 k = 2

Коэффициенты соответствующего ОДУ вычисляются через Я;- по формулам (10) или (12).

Формулы (і3) позволяют вычислить динамические характеристики колебательной системы (З)-декремент и частоту колебаний также через коэффициенты разностного уравнения (11).

4. Линейные дискретные модели динамических процессов и алгоритм вычисления коэффициентов ОДУ на их основе. В основе построения линейных дискретных моделей (ЛДМ) динамических процессов лежат известные решения ОДУ (3) и (З). Дискретизация непрерывных функций (4) и (7) с шагом t приводит к дискретным моделям вида:

X - X C^4Drk

xk =Xi %2Cer- , (1б)

1 - Ce4Dtk

xk = a0e -stk cos (otk + y0). (17)

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

Рассмотрим построение на основе выражений (1б) и (17) дискретных моделей, линейных относительно параметров, которые известным образом связаны с коэффициентами ОДУ. Из (1б)следует, что

, 4Dtk = Xk - Xi = ^ - X2Ce'f3t(k-i)

Ce

1 - Ce(k-i) '

Отсюда получаем хк-1 - хк-1 Хк—е _ х1 - х2 Х—ХL е .

Хк - Х2 Хк - Х2

После простых алгебраических преобразований следует

ХкХк-1 1 - е4Ъг )+ Хк (х1е"ГЪГ - Х2 )+ Хк-1 (х)2е4ЪТ - Х1 )_ Х1 Х2 (е^Т - 1).

Обозначив

Xk - X2

1б4

. _ ) - . _ Х - Х2е^Т • . (18.

11 _ і_^ 1 _ і_’ 13 _ ХіХ2’ (18)

получаем линейную дискретную модель, описывающую решение ОДУ (3):

1 Хк + Я2Хк-і +1 _ ХЛ-1- к _ 2,3,4... (19)

Нетрудно убедится, что параметры ОДУ вычисляются через коэффициенты ЛДМ (19) по формулам

П (2 - 1 )-У(11 + Я2 )2 + 413

(2 - 1 ) + V(( + І2 ) + 4Я3

a2 _ ■

t (( + Я2 ) + 4Яз

-, a1 _ -a2 (Я1 + Я2), a0 _ -a213 ■ (20)

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

xk = a0e_stk cos(atk + y0), xk_ = a0e-stVT cos(ank + y0 - at ) = a0e_stkesT [cos(fttk + y0 )cos at + sin(atk + y0)sin an], xk-2 = a0e_stke2st cos(ank + y0 _ 2wt ) = a0e~stke2st \cos(ank + y0 )cos 2wt + sin(ank + y0)sin 2an] с помощью простых алгебраических преобразований можно получить соотношение

1xk _1 _ Я xk _2 = xk, k = (21)

в котором коэффициенты связаны с параметрами ОДУ (6) по формулам

1 = 2e_sT cos an , 12 = e“2sT . (22)

Отсюда динамические характеристики системы (5) могут быть вычислены по формулам

1 1 о р ln 12

a = — arccos—-=, о =----------^=. (23)

n 2yj 12 arccos 2yj 12

Алгоритм вычисления коэффициентов ЛДМ (19) и (21) аналогичен рассмотренному выше на примере вычисления коэффициентов разностных уравнений (9) и (11).

Исследование помехозащищенности алгоритмов вычисления коэффициентов ОДУ при использовании разностных уравнений и линейных дискретных моделей. Целью исследований являлось изучение эффективности применения ЛДМ при вычислении коэффициен-

тов ОДУ. Исследования проводились с использованием персонального компьютера и программы Excel из Microsoft Office. При исследованиях с использованием математической модели, описываемой ОДУ первого порядка, рассматривалось дифференциальное уравнение

— = 2 _ 3x + x2 с начальными условиями x(0) = 3 . С различным шагом дискретизации t в со-

dt

/\ 4 _ et

ответствии с решением x(t) =----------- генерировалась выборка объемом N=10. Точные значения

2 _ e

округлялись до трех (е = 0,001) и двух (е = 0,01) знаков после запятой. Этим самым моделировалась помеха в результатах измерений при экспериментальных исследованиях. По выборке на основе разностного уравнения (9) и линейной дискретной модели (19) вычислялись коэффициенты ОДУ (3) и оценивалась относительная погрешность вычислений в %. По результатам исследований построены графики зависимостей погрешности вычисления коэффициентов ОДУ (рис.2- 4).

При исследованиях с использованием ОДУ второго порядка моделировалось решение уравнения x" + 0,1x’ + 4p2 x = 0 с начальными условиями, соответствующими a0 = 1. При этом начальная фаза y 0 изменялась случайным образом, а полученные результаты усреднялись по десяти различным значениям y 0 . Схема проведения исследований аналогична рассмотренной выше. По результатам исследований построены графики зависимостей погрешности вычисления динамических характеристик колебательной системы (5) (рис.5) от шага дискретизации при е = 0,001.

Проведенные исследования позволяют сделать следующие выводы.

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

0,02 0,03

0,04

0,05

0,06

Р и с. 2. Зависимости погрешности вычисления коэффициента а0 от периода дискретизации при использовании разностных схем (кривая 1) и линейных дискретных моделей (кривая 2)

Р и с. 4. Зависимости погрешности вычисления коэффициента а2 от периода дискретизации при использовании разностных схем (кривая 1) и линейных дискретных моделей (кривая 2)

0,02 0,03

0,04

0,05

0,06

Р и с. 3. Зависимости погрешности вычисления коэффициента а} от периода дискретизации при использовании разностных схем (кривая 1) и линейных дискретных моделей (кривая 2)

Р и с. 5. Зависимости погрешности вычисления декремента колебаний от периода дискретизации при использовании разностных схем (кривая 1) и линейных дискретных моделей (кривая 2)

т

т

тойчивостью вычислений при решении нормальной системы уравнений в МНК. Таким образом, при использовании разностных уравнений не удается получать результаты с высокой степенью точности.

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

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

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

БИБЛИОГРАФИЧЕСКИЙ СПИСОК

1. Агафонов С.А., Герман А.Д., Муратова Т.В. Дифференциальные уравнения: Учеб. для вузов/ Под ред.

В.С.Зарубина, А.П.Крищенко. М: Изд-во МГТУ им. Н.Э.Баумана, 1999. 336с.

2. Годунов С.К., Рябенький В.С. Разностные схемы. М.: Наука, 1973. 400с.

3. ЛинникЮ.В. Метод наименьших квадратов и основы теории обработки наблюдений. М.: Физматгиз, 1962.

320 с.