Диагностика и идентификация объектов в частотной области

Бондаренко Л.Н.

ДИАГНОСТИКА И ИДЕНТИФИКАЦИЯ ОБЪЕКТОВ В ЧАСТОТНОЙ ОБЛАСТИ

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

В ключевом докладе на 10-м международном симпозиуме Технического комитета ТС7 ИМЕКО [1] отмечено, что одним из основных применений измерений является построение и проверка достоверности моделей систем. Так как модели часто описываются значительным числом параметров, то для решения этих задач требуется повышение точности измерений, что может быть достигнуто в ближайшие годы за счет построения новых измерительных приборов на основе использования бурно развивающихся нанотехнологий.

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

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

Без ограничения общности можно считать, что в некотором диапазоне частот ] линейная двух-

полюсная электрическая схема замещения исследуемого объекта описывается функцией импеданса

(Р) = Ап(Р)!Вп(Р) , где Ап(р) = ^=0ап,кРк , В„(р) = '£1=0Ь„лрк , Ъпп = 1 , причем коэффициенты {апк

) п—1

неизвестны.

{Ь £.}П_о многочленов и структурный параметр

Пусть на различных частотах {тк^ , где N > п , ак > 0 , юк е[т1П,тах] измерены значения

2пк = %п(Рк) , где 2пк =ипк ^З^пк при рк =]тк . Тогда декомпозиция проблемы идентификации приводит к рассмотрению задач:

1) по измеренным на частотах {тк^ оценкам [ип , {Ул действительной и мнимой части значе-

ний функции (Р) требуется найти параметр П и оценки а , к }П=1о, Ь .*} П=1о коэффициентов {Эп к}к==о , {Ьп к }П=о при условиях выполнения заданной точности: | Эп к ~ Эп к I — еп к , \—пк _ —пк \—£п к ;

2) по найденным оценкам {ап^}п=о , {Ьпк}П=° найти оценки параметров элементов электрической схемы замещения.

При N > П первая из выделенных задач относится к аппроксимации рациональной функцией 7п(р) = Ап( Р)/ Вп (р) , р = ] т набора {тк . (иП.к. ^П,к )}^1 . Так как число неизвестных коэффициентов рациональной функции 2П меньше числа уравнений 2N , то естественным подходом к решению этой задачи является использование метода наименьших квадратов, непосредственное применение которого приводит к нелинейной системе уравнений. Поэтому в [3] для этого случая была предложена итерационная процедура, которая в дальнейшем многократно модифицировалась [4]. К существенным недостаткам этого подхода относятся значительная чувствительность к возмущению исходных данных и невозможность получения точного решения при N = П .

При N = П первая задача сводится к интерполяции рациональной функцией, а ее решение можно найти с помощью системы уравнений

Апатк)-

имеющей при определенных условиях единственное решение, не зависящее от набора испытательных частот {тк}" . Этот подход также называется методом конечно-частотной идентификации [4], а его существенным недостатком является плохая обусловленность матрицы системы (1). Поэтому для уменьшения числа обусловленности этой матрицы стараются соответствующим образом подобрать испытательные частоты

{т}П , что является достаточно сложной задачей. Следует отметить, что даже при полиномиальной интерполяции уже при числе уравнений больше четырех избегают напрямую решать соответствующую систему уравнений [5].

Рассмотрим иной подход к описанной задаче аппроксимации, базирующийся на методе рациональной интерполяции Кронекера [6] и сводящий ее решение к рекуррентной процедуре, в результате применения которой интерполяционный полином, построенный по исходным данным, последовательно трансформируется в

искомую ^ ( Р) .

Для этого при заданном частотном диапазоне [тт!п.ттах ] нормируем частоту т с помощью следующего преобразования 3 = й) / #?тах . Тогда набор испытательных частот заменится набором {Зк}^ , где

Зк <Е[у, 1] , У >0 и у = а)тт / й?тах . По исходным данным {Зк}^ и {^пк}к=^ можно построить интерполяционный полином Н 2д/ _1 (р), р = ]3 степени 2 N — 1 с вещественными коэффициентами, удовлетворяющий равенствам

ЯеН2Ы_,

Полином /“/ 2д/ (р) в соответствии с системой (2) принимает значения {1п к = ипк + }\/пк}к=-\ в интерполяционных узлах множества [—У. — У г] и [Уг, У] - состоящего из двух отрезков мнимой оси.

Рекуррентная процедура получения Zn (р) = Ап( р) / Вп(р) сводится к построению последовательности

многочленов соответствующих степеней, причем В-1 (Р) = 0 , В0(р) = 1,

Л2Л/-1 (Р) = Н 2Л/ -1 ( Р) - Л2Л/ (Р) = ^ 1\и~Р2 - Рк) ' где константу К можно, например, выбрать равной коэффициенту при старшей степени многочлена /~/ 2д/ _1 (р) . Эта процедура задается при к =0,1,..., N —2 соотношениями

Л2Ы-к-2(Р) = (акР + ^2Ы-к^(Р) ~ ^2Л/-/((Р)>

Вк+^р)={акр +рк)Вк{р)-Вк_^р).

На первом шаге ( к = 0 ) алгоритма Кронекера многочлен уА2д/ (р) в первом соотношении (3) делится на многочлен /\2д/_1 (р) . В результате чего находятся числа С^0, Д и многочлен ^2Л/-2 (Р) • Затем из второго соотношения (3) определяется в<|(р) . На следующих шагах алгоритма выполняется аналогичная последовательность действий.

Алгоритм Кронекера, по сути дела, основан на алгоритме деления многочлена на многочлен. При этом используется вычитание, которое во многих случаях приводит к сильному увеличению относительной ошибки. Это явление, например, особенно проявляется при использовании для построения многочленов

{^2Ы-к (Р)}о степенного базиса {Рк}о ' который является переполненным [5]. Применение этого базиса также является основной причиной плохой обусловленности систем уравнений (1) и (2).

Если в (3) выполняется неравенство 1 + с!едуА2д/ —/с—1 (Р) < /А2д/_к (р) , то описанный алгоритм Кронеке-

ра вырождается, так как невозможно определить параметры &к. 0 . В этом случае в (3) вместо двучлена

акр + рк следует ввести многочлен от рстепени >2 с неопределенными коэффициентами.

Для устранения недостатков степенного базиса описанный алгоритм Кронекера можно модифицировать за счет использования для записи многочленов {А2ы-к (Р))о базиса, построенного на основе многочленов

Чебышева - Ахиезера [7]. Считая в этом случае для простоты у достаточно малым, интерполяционный полином /~/ 2/\/(Р) г 1,1] надлежит разыскивать в следующей форме

" 2Л/-1 "

н2л/-1

к =0

где Ск(р) - многочлены Чебышева комплексного переменного, определяемые соотношениями: Ск(р) = сЬкв, р = У сИ<9, причем С0(р) =1 , С,(р) = -У р .

Свойства многочленов Чебышева позволяют по исходным данным находить коэффициенты {Ик }^ 1 в

по довольно простым формулам в двух важных случаях: {+]Зк}^ - узлы I рода, где

Зк = СОЗ((2Л/ - 2к + 1)лг / 4Л/) , и {±]Зк}^ - узлы II рода, где Зк = COS((N — к)7Г / (2Л/ — 1)) (в вещественном случае этот результат хорошо известен [5, 8], а эти системы узлов асимптотически оптимальны при

аппроксимации функции интерполяционным полиномом).

Отметим, что при вычислении коэффициентов полинома Н2д/_1 (р) по соответствующим формулам на узлах

I и II рода абсолютная ошибка не превосходит абсолютной ошибки исходных данных {и к^ , {V к}^ [5,

7] . Для рассматриваемых узлов величина у = 3^ , а узлы II рода имеют преимущество, заключающееся в

том, что для них Зы =1 . Узлы I рода являются нулями многочлена С2д/ (р) , а узлы II рода - нули многочлена С2д/ (р) — С2д/_2(р) * Поэтому при использовании в алгоритме Кронекера - Чебышева узлов I рода полагаем ^[2N (Р) ~ К

С2 N (Р), а в случае узлов II рода - ^2N (р) - К (С2Л/(Р) ^2Л/-г(Р)) •

Подчеркнем, что рекуррентная процедура Кронекера - Чебышева с вычислительной точки зрения не сложнее соответствующей процедуры Кронекера, так как при делении многочленов в базисе {Ск(р))о также можно использовать схему Горнера [8] . Это связано с тем, что многочлены Ск(р) = Тк(—] р) , а многочлены Чебышева Тк(3) удовлетворяют простому рекуррентному соотношению Т0(3) = 1 , Т^(3) = 3 ,

тк+2^&) = 2юГЛ+1(ю)-Тк{а>) , к = 0,1......

Описанный алгоритм Кронекера - Чебышева при N = П решает задачу рациональной интерполяции по исходным данным. Рациональная интерполяция обладает замечательным свойством (аналогичным свойству полиномиальной интерполяции): точная интерполяция по значениям заданной рациональной функции при числе интерполяционных узлов большем числа параметров этой функции дает в результате эту функцию. Это свойство можно использовать в алгоритме Кронекера - Чебышева при N > П для остановки работы процедуры.

В предположении нормировки частоты для нахождения по исходным данным {U n,*}N , К, * }N коэффициентов {^ }^ 1 интерполяционного многочлена (4) на узлах I и II рода в пакете аналитических вычислений

Maple была разработана соответствующая процедура. Также в этом пакете написана процедура, реализующая метод Кронекера - Чебышева в вещественной и комплексной областях на узлах I и II рода. В этой процедуре предусмотрено задание параметра идентификации tol (от слова tolerable - допустимый), связанного с величинами | Эпк — Зп)< \ — sn к , подбор которого позволяет определять структурный параметр n при фиксированном значении N > П . Также в этих процедурах была учтена возможность измерения значения функции импеданса на постоянном токе. В этом случае к набору испытательных частот К)" добавляется частота CDq = 0 . В процедурах это достигается, если за узлы I рода взять нули многочлена C2NAP) ' а за Узлы И рода - нули многочлена C2/v+i(P) — ^2/v-l(P) •

Отметим, что на практике испытательные частоты {<3^ могут устанавливаться только с определенной

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

Разработанные процедуры позволили на основе имитационного моделирования в пакете Maple протестировать большое число задач. В частности, были рассмотрены примеры аппроксимации методом Кронекера -Чебышева объектов со следующими функциями импеданса:

ч 1 ,0016р4 +15,0201 р3 + 77,0854р2 +167,1461 р +140,0858

а) Z( р) =---------—г—!-----^^-------------------’-----------’----.

р5 +17 р4 +107 р3 + 321 р2 + 474 р + 280

Разложение этой функции импеданса на простейшие дроби имеет вид

ч 1 3• 10—4 5• 10—4 8• 10—4

Z (р) =--- +-------Т" +---Z~ +----2—л---ТТ '

р + 2 р + 4 р + 7 р2 + 4 р + 5

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

Z 1( р) = р+2 .

Возможный диапазон изменения параметра идентификации в рассматриваемом случае оказался равным tol е [0,07; 0,6] . Для различного числа N испытательных нормированных частот I рода были получены следующие результаты при вычислении с шестью значащими цифрами:

" 1 00126 " 1,00087

N = 1 , Z 1( р) = ' ; N = 2 , Z 1( р) =  ,

р + 2,00135 1 р + 2,00047

о 1,00047 " 1,00655

N = 3 , Z 1( р) =------------ ; N = 4 , Z1( р) =--;

' 1 р +1,99957 ' 1 р + 2,01316

а, с \ 1,16554

N = 5 , Z 1( р) =—--------;

1 р + 2,36842

при вычислении с десятью значащими цифрами:

л, „ -,",4 1,000440343 „ ‘ ч 1,000449161

N = 4 , Z,

р +1,999505947 1 р +1,999525650

При выборе параметра идентификации № = 0, N = 1 , &$ = 0 и вычислении с шестью значащими цифрами

найдено

' , 0,919794 • 10~4 р +1,00109

2 Л р) = ---------------1----.

1 р + 2,00095

ч 3, 0013 р4 + 41 ,0162 р3 +187, 0689 р2 + 369,1158 р + 280, 0648

б) 2 (р) = ------—^^--------------------------’-------------------’-----.

р5 +17 р4 +107 р3 + 321 р2 + 474 р + 280

Разложение этой функции импеданса на простейшие дроби имеет вид

, 1 2 5• 10~4 8• 10~4

2( р) =-----П +-----л +---т" +-2----------"л-сГ ,

р + 2 р + 4 р + 7 р2 + 4 р + 5

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

3 р + 8 = р2 +6р + 8 '

Возможный диапазон изменения параметра идентификации в рассматриваемом случае оказался равным № € [0,0003; 0,09] . Для различного числа N испытательных нормированных частот I рода были получены следующие результаты при вычислении с шестью значащими цифрами:

2 68770  3,00106 р + 7,99359

N = 1 , 22(р) = 2,68770 ; N = 2, 22(р) =-. н -

р + 2,70667 2 р2 + 5, 99641 р + 7,99170

3, 02693 р + 8, 30287

N = 3 , Z 2( р) = _ ;

р2 + 6,13776 р + 8, 29835

при вычислении с десятью значащими цифрами:

N

N

3 ,

4 ,

Z 2( p) =

Z 2( p) =

N = 5 , Z2( p) =

3, 0001152653 p + 7, 984270576 p2 + 5, 992020527 p + 7, 982509894 3, 001667312 p + 7, 998034916 p2 + 5, 998688843 p + 7, 996088053 2, 989658263 p + 7, 895311448

р2 + 5, 948187750 р + 7, 894884429

При выборе параметра идентификации tol = 0, N = 1 ,

0 и вычислении с шестью значащими цифрами

—0, 729146 • 10—3 p2 + 3, 00575 p + 8, 07257

Z 2( p) =

р2 + 6,03068 р + 8,07077

при выборе параметра идентификации tol = 0, N = 1 ,

ми найдено

: 0 и вычислении с десятью значащими цифра-

0, 7227776518 • 10—4 р2 + 3, 000632674 p + 7, 989725939

Z2( p) =

р2 + 5, 994524628 р + 7, 987877316

При задании параметра идентификации tol = 0, N = 5 и вычислении с двадцатью значащими цифрами заданные функции импеданса были восстановлены с большой точностью. Аналогичные результаты также получены для испытательных частот II рода.

Достоинствами предлагаемого метода аппроксимации Кронекера - Чебышева является как простота вычислительных процедур, так и уменьшение чувствительности к ошибкам исходных данных, причем последнее

можно усилить за счет разумного выбора частотного диапазона [^mjn, ^max] .

Так как не найдены аналоги этого метода в комплексной области, то для сравнения рассматривались примеры рациональной аппроксимации в вещественной области тригонометрических и гиперболических функций на отрезке [tu / 2], показательной функции на отрезке [0,1] , функций Бесселя и т. п.

Для этого использовались разработанные процедуры и в высшей степени оптимизированные стандартные процедуры пакета Maple, реализующие метод Ремеза наилучшей аппроксимации и метод аппроксимации Паде - Чебышева, разработанный Кленшоу и Лордом [6].

На рис.1 показаны графики логарифмической погрешности рациональной аппроксимации функции Бесселя

Jq( X) на отрезке [0, 7о] , где j * 2 40483 - наименьший положительный нуль функции Jq(x), с использованием алгоритмов Ремеза, Паде - Чебышева и Кронекера - Чебышева на узлах I и II рода.

Рис.1. Погрешность аппроксимации функции Бесселя J0(X) .

На этом рисунке Igerror = lg max | J0(x) -Rn(x) |, где аппроксимирующая рацион

Xе[0,j0]

альная

функция Rn

имеет соответственно степень числителя [П /2] и степень знаменателя [(n +1) / 2] , а [а] означает целую часть числа а . По рис. 1 видно, что методы Ремеза, Паде - Чебышева и Кронекера - Чебышева по качеству аппроксимации практически равноценны, но метод Ремеза работает по времени в десятки раз дольше, чем остальные рассматриваемые методы.

Таким образом, ближайшим конкурентом метода Кронекера - Чебышева в вещественной области является метод Паде - Чебышева, но для его работы необходим отрезок ряда Фурье - Чебышева аппроксимируемой функции.

Отметим, что вторая задача, полученная в результате декомпозиции проблемы идентификации, обычно называется задачей реализации и проще задачи аппроксимации. Эффективные методы ее решения основаны на применении канонических форм электрических двухполюсников [9]. На базе этих методов в пакете Maple разработана универсальная программа для решения задачи реализации, использующая известные канонические формы Фостера, Кауэра и Ли, а также две новые канонические формы.

ЛИТЕРАТУРА

1. Финкелстайн Л. Измерения в мягких системах // Датчики и Системы. - 2004. - № 10. - С. 61-67.

2. Агамалов Ю. Р., Бобылев Д. А., Кнеллер В. Ю. Виртуальные измерители-анализаторы параметров импеданса // Датчики и Системы. - 2004. - № 5. - С. 14-18.

3. Sanathanan C. K., Koerner J. Transfer function synthesis as a ratio of two complex polynomials // IEEE Trans. Automat. Control. - 1963. - V. AC-9. - № 1. - P. 56-58.

найдено

4. Орлов Ю. Ф. Идентификация по частотным параметрам при параллельных испытаниях // Автоматика и телемеханика. - 2007. - Т. 68. - № 1. - С. 20-40.

5. Бабенко К. И. Основы численного анализа. - М.: Наука, 1986.

6. Бейкер Дж., Грейвс-Моррис П. Аппроксимации Паде. - М.: Мир, 1986.

7. Бондаренко Л. Н. Определение параметров передаточной функции средств измерений по значениям амплитудно-частотной и фазочастотной характеристик // Датчики и системы. - 2004. - № 7. - С. 18-20.

8. Люк Ю. Специальные математические функции и их аппроксимации. - М.: Мир, 1980.

9. Бондаренко Л. Н. Структурно-параметрическая идентификация многоэлементных двухполюсников //

Информационная измерительная техника: Межвуз. сб. науч. тр. - Вып. 29. - Пенза: Информационно-

издательский центр, 2005. - С. 40-45.