Помехозащищенный метод определения параметров линейной динамической системы на основе импульсной характеристики Текст научной статьи по специальности «Математика»

ВЕСТН. САМАР. ГОС. ТЕХН. УН-ТА. СЕР. ФИЗИКО-МАТЕМАТИЧЕСКИЕ НАУКИ. 2007. № 1(14)

Математическое моделирование

УДК 519.246 В. Е. Зотеев

ПОМЕХОЗАЩИЩЕННЫЙ МЕТОД ОПРЕДЕЛЕНИЯ ПАРАМЕТРОВ ЛИНЕЙНОЙ ДИНАМИЧЕСКОЙ СИСТЕМЫ НА ОСНОВЕ ИМПУЛЬСНОЙ ХАРАКТЕРИСТИКИ

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

Проблема параметрической идентификации линейной динамической системы является одной из важнейших в математическом моделировании. При описании линейного динамического объекта в задачах управления, оптимизации или прогнозирования обычно используется его передаточная функция W (p), которая, по сути, представляет собой изображение импульсной переходной (весовой) функции — реакции системы на импульсное (ударное) воздействие [1, 2]. Динамические характеристики линейной системы, непосредственно связанные с коэффициентами соответствующего однородного дифференциального уравнения, вполне могут быть описаны через параметры импульсной переходной функции (импульсной характеристики системы). Поэтому среди известных методов параметрической идентификации особое место занимают относительно простые, но, в то же время, эффективные методы, в основе которых лежит измерение мгновенных значений импульсной характеристики системы. В частности, при определении динамических характеристик колебательных систем широко применяется метод затухающих колебаний, который заключается в записи виброграмм свободных колебаний диссипа-тивной системы и определении по ним логарифмического декремента колебаний [3]. Однако часто алгоритмы, используемые при определении параметров линейной динамической системы по мгновенным значениям импульсной или переходной характеристики, не ориентированы на применение статистических методов обработки результатов эксперимента. Вследствие этого из-за случайной помехи в результатах наблюдений эти методы обладают низкой помехозащищенностью [4]. Решить проблему качественного изменения способов определения параметров линейной системы по результатам измерения мгновенных значений ее импульсной характеристики можно на основе линейно параметрических дискретных моделей в форме стохастических разностных уравнений [5].

Пусть динамический процесс ~ (t) в линейной системе при входном воздействии x(t) описывается дифференциальным уравнением с действительными коэффициентами вида ^ dj у (t) djx(t) ^

у cj-— = у bj-—. При импульсном воздействии на систему, которое описывается

j=o dt j=o dt

8 -функцией Дирака, реакция системы — весовая функция (импульсная характеристика) — при отсутствии кратных корней соответствующего характеристического уравнения может быть

р

представлена в виде у(() = у a¡ exp(a¡t), где в общем случае a¡ = Rea,. + i Im a¡ и

i=1

a¡ = Re a¡ + i Im a¡ — комплексные числа. Так как коэффициенты дифференциального уравнения действительные числа, то для каждого комплексного корня характеристического уравнения ai найдется сопряженный ему корень a¡, причем соответствующие им коэффициенты a¡ и a¡ также являются комплексно-сопряженными числами. В таких случаях сумма a¡ exp (a¡t)+ a¡ exp (a¡t) описывает затухающие колебания и может быть представлена в виде

Д. exp (Re att) cos (im att ), где начальные амплитуда At = 2^/Re2 at + Im2 at и фаза Im a..

У =

arctg

Re a

Re a . > 0,

Im a

p + arctg-Re at < 0

Re a

— действительные числа. По найденным оценкам динамических

характеристик а1 коэффициенты дифференциального уравнения с., . = 0,1,..., р -1, с = 1,

могут быть определены по формулам: ер-1 = -Z

c=(-1)p-1 Zaa• a<Pi' co =(-1))П a •

a i' cp-2 = Za

i, j=1

aj . CP-3 = -Z i, j,k

i<j<k

a.a a k

>1. '2.' >p-1 I1<I,<- <iP_

Рассмотрим задачу оценки параметров а, и а, (,' = 1,2,..., р ) функции ~(() = Хя, ехр(а/)

/=1

по результатам измерений ук = ~(хк) + ек ее мгновенных значений, где % - период дискретизации; ек — случайная помеха в результатах наблюдений; к = 0,1,2,* , N — 1; N — объем выборки результатов измерений.

Задача решается на основе линейно параметрической дискретной модели, связывающей р +1 последовательное значение дискретной функции

у

ук=Х я ехр(а%к К

/=1

р

которую можно представить в виде ук = X а/' , где //., = ехр (а,%).

(1)

Применяя к дискретной функции (1) г -преобразование, имеем 2{ук }=Х а1-—. От-

,=1 1 — /,г

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

Лр(х-1)2{ук } = Вр—, (2)

k p -p-1

где Лр(г 1 )=П( — /¡г 1 )= 1 — X и Вр—1 ( 1 )=ХЬ,г ' — многочлены степени р и р — 1 от-

,=1 '=1 '=0

носительно переменной г -1 . Очевидно, что коэффициенты многочлена Лр (г-1 ) описываются

формулами

У У

1 =Z m>. 1 = -Z mm j.

i=1 i,j

i< j

p

яр-1 =(- 1)P Z m, mh • mip-1.

У

1 = Z mm j mk. •

i, j ,k i< j<k

1P = (-1))+1 П m. •

(з)

>l h" ■ ip-1

£< i, <• < ip -

Из (2) с помощью обратного г - преобразования можно получить

~к — КУы—1 — ^к—2 — * — Яр~к—р = К8к + ЬА—1 + Мк—2 + * + Ьр—18к—(р—1) ,

1, при к = 0,

(4)

где 8к = < — символ Кронекера, к = 0,1,..., N — 1.

[0, при к Ф 0

С учетом (1) введем обозначения

р р р

яр+1 = й =X а, яР+2 = й = X ат<, 1р+3 = й = X ат,2

яр+j = ~j-1 =Z a mi-1>

1 p = ~p-1 = Z am"-1.

(5)

i=1

¡=1

i=1

i=1

i=1

i=1

i=1

Используя (4) и (5), коэффициенты многочлена Bp-1 (z 1) можно представить в виде

Ь0 = Яр+1 , Ь1 = Яр+2 - АЯр+1 , Ь2 = Яр+3 - Я1Яр+2 - Я2Яр+1 , ' ,

1

Ъ1 = Яр+1+1 - ЛЯр+1 - Я2Яр+1-1 — - ЯА+1 = Яр+1+1 ЯА+1+1-1, * , (6)

1=1

Р-1

Ьр-1 = Я2р - Я1Я2р-1 - Я2Я2р-2 - * - Яр-1Яр+1 = Я2р - X Я1Я2р-1 •

1=1

При этом линейно параметрическая дискретная модель, связывающая в виде рекуррентной формулы р +1 мгновенное значение дискретной функции (1), принимает вид

к = 0,1,..., р -1,

1 У к-1 + 12 У к-2 + * + 1 рУк - р , к = Р, р +1,•••, N -1,

(7)

где коэффициенты 1 ^ (] = 1, 2,2р), определяются соотношениями (3) и (5).

При экспериментальной обработке импульсной характеристики формируется выборка результатов измерений ук = ~к + ек, которые содержат случайную помеху ек, к = 0,1,..., N -1, N — объем выборки. В этом случае линейно параметрическая дискретная модель (7) принимает вид стохастического разностного уравнения

\+1+к +ек, к = 0,1,...> р -1 р 1 (8)

X }^Ук -1- X1 е к -1+ек, к=p, р+и.,N -1 1=1 1=р

Ук =

или в форме обобщенной регрессионной модели:

Гй = F1 + h,

р (9)

Здесь 1 = (11,12, * , 12 )т — вектор неизвестных коэффициентов линейно параметрической дискретной модели; е = (е0,е1,...,еN-1)T — N-мерный вектор случайной помехи в результатах наблюдений; " = (Пр"П^—,)Т — N -мерный вектор эквивалентного случайного возмущения в стохастическом разностном уравнении; Ь = (у0, у1,..., УN-1)T — N -мерный вектор правой части; F = [/1 • /2 •. • /' • / •. • /2 ] — матрица регрессоров размера N х 2р , столбцы которой описываются формулами: / = (°Д» ,0, Уp-l, Ур, УN-2)T, ./2 = (0А* ,0, Уp-2, УР-х' , УN-з)T, * ,

/ = (0,0,- ,0, У0, У1,* , УN - p-l)т, /р+1 = (1,0,- ,0)\ /р+2 = (0,1,- ,0)\ • , /2 р = (0,0,- ,1,0,- ,0)\

Матрица Р размера N х N в стохастическом разностном уравнении эквивалентного возмущения — нижняя треугольная. Первые р строк матрицы: р1 =(1,0,0,* ,0),

р2 =(0,1,0,* ,0), * , рр =(0,0,* ,1,0,- ,0) — описываются формулой р9 =Ц, 1 '

1 = 1,2,..., р,, 1 = 1,2,..., N. Остальные строки матрицы Р имеют вид:

Р,

= (-1p,-lp-1,-lp-2,' , "11,1,0,. ,0) , Pp+2 =(0,-1p, -1p-„-lp-2,' ,-12,-11,1,0,« ,0),

p+1 \ p> p-1> p-2> > 2> j'-'/? rp+2 Y"'? /up' /up-1> p-2 >

pw =(0,. ,0,-1p,-lp-1,-1p-2,. ,-12,-1i,1).

В предположении, что случайные возмущения eк имеют нулевое математическое ожидание, одинаковые дисперсии и не коррелированны между собой, помехоустойчивое определение параметров импульсной характеристики ai и ai сводится к среднеквадратичному оцениванию коэффициентов 1 линейно параметрической дискретной модели (8).

При среднеквадратичном оценивании коэффициентов 1j могут быть использованы различные схемы вычислений. Наиболее простым является алгоритм, в основе которого лежит 140

I |2 N ||2

минимизация функционала h = b - Fk\\ ^ min. МНК-оценки коэффициентов стохастического разностного уравнения в этом случае вычисляются по формуле

X =(FTF)-1 FTb. (10)

Однако такой подход мало эффективен из-за большого смещения оценок, обусловленного корреляцией между элементами hk эквивалентного случайного возмущения. Другой, более эффективный и помехозащищенный, метод среднеквадратичного оценивания коэффициентов стохастического разностного уравнения основан на минимизации функционала ||e|2 =|\y-|2 ^min. Главной особенностью такого подхода является применение итерационной процедуры уточнения среднеквадратичных оценок коэффициентов Xj. Это позволяет практически устранить

смещение в оценках и, тем самым, добиться высокой точности вычисления параметров исследуемой функции [5].

Рассмотрим алгоритм итерационного метода среднеквадратичного оценивания коэффициентов разностного уравнения. Из формулы (9) получаем P-lb = P— Fl + e . Элементы p- обрат-

1, j= , 0, j * i.

ной матрицы РЛ1 зависят от коэффициентов % и описываются формулами p- =

при i = 1,2,..., p , j = 1,2,..., N и p- =

p

1+E Wk, j, j=и

k=1 p

E ^p-k,j, j *i

k=1

при i = p +1, p + 2,..., N,

j = 1, 2,..., N.

В основе итерационной процедуры лежит минимизация функционала ||e||2 » ¡Р.-) - Px-k1)Fl|| ^ min , где P-' - обратная матрица, при вычислении элементов которой используются оценки Xj) коэффициентов разностного уравнения, найденные на k -ой итерации. Очевидно, что данный функционал представляет собой квадратичную форму относительно искомых коэффициентов Xj. Следовательно, он достигает своего неотрицательного минимума. При этом нетрудно показать, что минимум данного функционала достигается в точке

€(k+1)=(FT Q^F)-1 FTQ^k) b, (11)

где Q,4=( P-1)1 (P-1).

Итерационная процедура среднеквадратичного оценивания коэффициентов разностного уравнения состоит в следующем. На первом шаге по формуле (10) находится начальное приближение €(0) = (FT F)"' FTb . Затем, полагая в формуле (11) k = 0 и %к' = Х0, вычисляется следующее приближение €(1)=(FтQ^ F) FтQ^ b. Оно вновь подставляется в правую часть

формулы (11) и находится новое приближение Х2" и т.д. Процесс уточнения среднеквадратичных оценок коэффициентов стохастического разностного уравнения продолжается до тех пор,

пока выполняется условие ||Xk+1) - Xk0,01||Xk. Результаты численно-аналитических исследований показали хорошую сходимость итерационной процедуры (не более десяти итераций). Найденные среднеквадратические оценки Xj коэффициентов линейно параметрической

дискретной модели (8) лежат в основе вычисления помехозащищенных оценок параметров импульсной характеристики ai и ai, i = 1, 2,..., p .

Для этого сначала находятся корни im, i = 1, 2,..., p , алгебраического уравнения

mp + f1mp-1 + f2mp-2 + • +fp-1m+fp = 0. (12)

Затем по формулам d€ = t l\nвычисляются оценки параметров ai (i = 1, 2,..., p ). На заключительном этапе по найденным и im посредством решения системы линейных алгебраиче-

141

ских уравнении

к + к + • +а ь = 4+1, 1+2,

М +¡12&2 +• + ! РкР =

¡2а1 +¡22 а + +¡4 = = 4+3,

¡Р-1 к + ¡р -1а +• + ¡р- = р

вычисляются оценки коэффициентов а1 (1 = 1, 2,..., р ) в функции (1).

При комплексно-сопряженных корнях алгебраического уравнения (12) ¡¡к = Яе ¡¡к +' 1т ¡¡к и ¡к+1 = ¡к = Яе ¡¡к -11т ¡¡к и соответствующих им комплексно-сопряженных коэффициентах ак = Яе к +11т к и к+1 = к = Яе к -11т к, наИденных из системы уравнении (13), пара слагаемых к ехр((кк/) + к+1 ехр (акк+1/) в математической модели импульсной характеристики р

) = ^ к ехр (кк1), описывающая затухающие колебания, может быть представлена в виде

А ехр (Яе (кк/) 008 (1т скк/+укк),

(14)

где Яе&к = ' х Яе2 ¡1к + 1т2 ¡1к и = 1т(£к =

arotg

р + arotg

1т тк Яе тк' 1т тк

Яе тк

Яе тк > о,

Яе т к < о

arotg

— оценки показателя

затухания и частоты колебаний, Д = 2д/Яе2 к + 1т2

р + arotg

1т к

Як"'

1т ак

Яе ак

Яе к > о, , Яе а < 0

оценки начальной амплитуды и фазы колебаний.

Таким образом, в работе получены следующие результаты.

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

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

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

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

и

1. Эйкхофф П., Ванечек А., Савараги Е. и др. Современные методы идентификации систем / Под ред. П. Эйкхоффа. М.: Мир, 1983. 400 с.

2. Штейнберг Ш. Е. Идентификация в системах управления. М.: Энергоатомиздат, 1987. 80 с.

3. Писаренко Г. С., Яковлев А. П., Матвеев В. В. Вибропоглощающие свойства конструкционных материалов: Справочник. Киев: Наук. думка, 1971. 376 с.

4. Мелентьев В. С. Методы и средства измерения параметров электрических цепей на постоянном токе. Самара: СамГТУ, 2004. 120 с.

5. Зотеев В. Е., Попова Д. Н. Определение динамических характеристик нелинейных диссипативных систем на основе стохастического разностного уравнения // Вестн. Сам. гос. техн. ун-та. Сер.: Физ. - мат. науки, 2006. Вып. 42. С. 162-168.

6. Семёнычев В. К., Зотеев В. Е. Определение параметров затухающих колебаний на основе разностных схем // Проблемы прочности, 1988. № 12. С. 101-105.

Поступила 15.07.2006.