Идентификация импеданса Варбурга на базе RC-схем замещения

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

ИДЕНТИФИКАЦИЯ ИМПЕДАНСА ВАРБУРГА НА БАЗЕ RC - СХЕМ ЗАМЕЩЕНИЯ

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

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

С этих позиций построение электрической схемы замещения некоторого объекта по результатам измерений значений его амплитудно-частотной характеристики (АЧХ) и фазочастотной характеристики

(ФЧХ) на фиксированных частотах {Фу}—^ в заданном диапазоне ,Фтах] является задачей струк-

турно-параметрической идентификации в частотной области.

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

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

Так как идентификация представляет собой итерационный процесс построения модели объекта, то использование линейной модели на первом шаге позволяет существенно упростить задачу структурнопараметрической идентификации без уменьшения ее общности.

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

Таким примером может служить задача построения модели импеданса Варбурга в виде конечного RC-двухполюсника. Этот импеданс является базовым элементом в схемах замещения, применяемых в методе переменноточной полярографии. Он описывает систему с распределенными параметрами, моделируемую импедансом полубесконечного ЯС-кабеля (трансмиссионной линии). Импеданс Варбурга используется в кинетике электродных процессов при анализе стадии разряда - ионизации и определяется следующими выражениями:

ZD = (1 - Wa~1/2 = RD - KCD)~1, ZD =42Wp~V1, p = jrn,

—1/2 —1 —1/2

где диффузионное сопротивление R^ =W-Ф и диффузионная емкость = W -Ф связаны

равенством RjjCjj Ф = 1 , а ю - круговая частота [2].

—1/2

Вместо импеданса Варбурга для простоты будем рассматривать функцию Z(Р) = р и- 2 , значения которой Z(py) = и у + jVy на практике находятся путем измерений на различных частотах {y}N ,

i i ЮттФтах], Фу У 0 . Все {(u у v)}^ являются вещественными, а время измерения с применением стандартных испытательных сигналов не превышает время переходных процессов.

Пусть Z(p) = Z(p) , P = j& получена нормировкой частоты (О = (О / £%1ах ' где е1Т; 1] i-У = фтт / фтах . Тогда поставленная задача заключается в определении по результатам измерений коэффициентов рациональной функции

_ Ап-\(р) „ п-1  _ п _  _

zn

t>n\P) i=0 ’ ¿=0 ’

аппроксимирующей при N>n функцию Z(/?) = Z(/>) и реализуемой в виде RC- двухполюсника с 2п элементами, причем следует подчеркнуть, что структурный параметр n также подлежит определению.

f ~ -*N

Применение относительных частот {Фул позволяет упростить обработку результатов, а фиксированные частоты {ф}^ используются только при измерениях. Так как измеряемые величины {(u?-, Vy )}N могут быть искажены за счет влияния внешних факторов, а значения функции Z СР) легко вычисляются

,p N

в точках {Фу} , то с помощью имитационного моделирования можно анализировать только погрешности,

получаемые при обработке "результатов измерений".

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

—1/2

Покажем, что функцию f (z) = z можно с любой степенью точности аппроксимировать аналитиче-

ской функцией F (z) класса S , регулярной во всей комплексной плоскости с разрезом вдоль полуоси (—&, 0) , Im F (z)/Im z < 0 при Im z Ф 0 и допускающей представление [4]

“ da

F

0 0 z +t

где cr(t) (0 < t <rc) - некоторая неубывающая функция, J(0(1 +1) 1d^(t) < 0 , а Pq — 0 . Это вытекает из

следующей формулы преобразования Стилтьеса [5]

*z—1/2e«zerfc(a1/2z1/2) = ?t~Vle~atdt , 0 t + z

где Re^>0 , erfc(w) = 1 — erf(w) , а erf(w) - функция ошибок.

В [4] доказано, что (2) аппроксимируется рациональным выражением

п Pi

Sn

и i=\z + ^ 1 1

Реализация этого выражения осуществляется ЯС-двухполюсником [3].

Поэтому естественно использовать задачу Золотарева Е. И. о нахождении такой рациональной функ-1/2

ции (3), что величина z Sn (z) наименее уклоняется от единицы на данном вещественном отрезке —2

[1,к 2] . При Р0 = 0 > 0 решение этой задачи дается в эллиптических функциях соотношениями [6]:

n—1

0 ч 2M1 Д

Sn

1 + 1 j

i=1 1

, л2,„, , ,,П 1 - к' 2 sn2(K' / n; к )sn2((2i - 1)K' / n; к) _

1 = dn

i=1 1 — к 2 sn2

где sn(w;к), cn(w;к), dn(w;к) - эллиптические функции Якоби модуля к, к’ = >/1 — к2 - дополнительный

модуль, а K - полный эллиптический интеграл первого рода для модуля к .

Ограничения р0=° ^ >0 приняты для определенности, чтобы число параметров идентифицируемого

RC-двухполюсника было равно 2n , но все соотношения могут быть модифицированы для нечетного числа параметров.

Применяя линейную замену переменного z, можно легко перейти к требуемому действительному от-

2 —1/2 резку [У,1], где у = к . Тогда относительная погрешность аппроксимации функции f (z) = z при

фиксированном n с помощью выражений (4) - (5) на отрезке [У,1] будет минимальна и равна

S = max |1 — z1/2Sn

ze[y,1]

s, N —1/2 tv ч —1/2

Переход от функции f (z) = z к Z(p) = p осуществляется заменой z на p = ]Ф . Полагая для

- 1/2 ~

простоты — У-> ^тах = 1 f получим Z(p) = p , a p<E.\jy^j\. В этом случае приходится сравнивать

не только

АЧХ Щр)\ и ее аппроксимацию \znm. но и ФЧХ Ф() = arctg(Im Z(p) / Re Z(p)) и ее аппроксимацию Фп(3) = arctg(Im Zn(p)/Re Zn(p)) , где Zn(p) = Sn(p) .

Для этого используются следующие относительные погрешности

*А = шах "^1-1 ^)», ¿ф = max "V" . (7,

А й)е[у,1] \Z(p)\ в)ф,1] |ф(©)|

В большинстве случаев задача рациональной аппроксимации действительной функции на заданном отрезке решается приближенными методами. Довольно эффективным подходом к рациональной аппроксимации является метод аппроксимаций Паде-Чебышева (ПЧ). Для формального разложения функции f (z) на отрезке по многочленам Чебышева T-(z) = cos i arccos z

x>

f (z) = £ diTi (z) i=1

задача об ПЧ аппроксимациях состоит в построении рациональной функции, степени числителя и знаменателя которой равны соответственно m и п, с использованием первых n + m +1 коэффициентов разложения.

Удачный с вычислительной точки зрения алгоритм такой аппроксимации разработан Кленшоу и Лордом в 1974 г. [7], а его модификация реализована в пакете Maple в виде процедуры "chebpade". С помо—1/2

щью этой процедуры несложно найти рациональную аппроксимацию функции f (z) = z на отрезке [у,1],

—1/2 - г . .п

а также аппроксимацию функции ИКр) = р при pE-yjy^j].

V/-4 —1/2

Существенными недостатками рассмотренных подходов к аппроксимации функции ИКр)=р являются два момента: отсутствие связи рациональной аппроксимации с результатами измерений {(^, Vi )}^ на

частотах {^ }^ и невозможность определения структурного параметра п.

Поэтому для решения задачи структурно-параметрической идентификации был предложен алгоритм

Кронекера-Чебышева (КЧ) [8-10] . В алгоритме КЧ по данным {щ}\ и {(му>v/)}^ строится интерполяционный полином

2N-1

H1N_x{p)= X hftip), (8) i=0

где {CytJ?)}^ , /?е[— j,j] ~ многочлены Чебышева комплексного переменного, определяемые парамет-

рическими соотношениями С-{р) = COS/'# , р = jcosO .

Полином (8) принимает 2N значений {uy± jv^}N в соответствующих интерполяционных узлах {+jcbj}^ множества [— у,— jy\ UL/'/s./] г состоящего из двух отрезков мнимой оси, а его значения в этих

узлах удовлетворяют системе уравнений: = э 1^/^ = vi5

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

а >2N—1 ,0,

{hy }q в выражении (8) в двух частных случаях, один из которых характеризуется тем, что интерполяционные узлы являются нулями многочлена C^j^jip) г Т. е. Фу = COS((2A^ -2i + l)7Ul 4А0 •

В этом случае справедливо неравенство у< COS((2N — 1)W(4N)) , а испытательные частоты легко находятся a priori по формуле Фу = •

Рекуррентная процедура построения функции (1) в алгоритме КЧ состоит в вычислении последовательности многочленов {^2ЛГ—/0^ в ^азисе и послеД°ва1Т1ельнос1Т1И в степенном

базисе {р1} о- где в_^(р) = о, 2?q(J?) = 1, ^2TV-1^ =H2N-\(P) ' A2N^= C2N^ ' a также вспомога-

тельных параметров {(Ху, py}N 1 по формулам

A2N-i-2^ = {ai~P + Pi)A2N-i-№ - A2N-i

На первом шаге алгоритма КЧ полагаем z'=0 и делим многочлен в равенстве (9) на

* -^Ри Э1ТОМ находятся коэффициенты CCq, /3q и многочлен ^2N—2' а за1тем с пом°Щью выражения (10) определяется многочлен В^(р) . На последующих шагах алгоритма при / = 1,2,...,N — 1 выпол-

няется аналогичная последовательность действий.

При у> cos((2N — 1)W(4N)) алгоритм КЧ заменяется алгоритмом Кронекера-Чебышева-Ахиезера (КЧА) [8-10]. В алгоритме КЧА вместо многочленов Чебышева Cj(p) применяются модифицированные многочлены Чебышева-Ахиезера CAj(p) , задаваемые при р е [—j,— jy] ^\jy^j\ соотношениями

СЛ2m^ = Cm(i2p2+l + r2)l(}-y2)) , СА2т+\^ = РСЛ2т'Р) ’ ОТ = 0Л>----

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

В алгоритме КЧА — СА 2дК-Р) ' а ^2N—1^) находится по исходным данным {фу}^ и

{(ui,vi )}iN в виде полинома H2N-l(P) по системе модифицированных многочленов Чебышева-Ахиезера

{С4,-(р)}£

2N-1

h2N-i(p)= X hfAi(P)- (И) i=0

Построение функции (1) в алгоритме КЧА производится аналогично алгоритму КЧ, но для получения последовательности многочленов {А-2^_у(р)}^+^ г и нахождения коэффициентов {//у}2^ ^ поли-

нома (11) по отдельности рассматриваются четные и нечетные части многочленов в формулах (9) -(11).

В алгоритме КЧА интерполяционные узлы являются нулями модифицированного многочлена

Чебышева-Ахиезера г а испытательные частоты со- = Фшах^/' находятся a priori с помощью со-

отношения

- Il-y2 f (2N - 2/ + 1)яЛ 1 + r2

a>j=A----—cos --------------------------------------------— + — .

' 2 { 2 N J 2

В этой формуле нули G)j вырождаются в нули =C0S((2iV — 2i + Х)7Г / AN) обычных многочленов Чебыше-

ва N^ip) при у = 0.

Основой алгоритмов КЧ и КЧА является деление многочлена на многочлен. Так как многочлены Q(_p) и СА-{р) удовлетворяют простым рекуррентным соотношениям, то при его выполнении можно применять

вычислительную схему Горнера. Это использовано при реализации алгоритмов КЧ и КЧА в виде простых вычислительных процедур в пакете аналитических вычислений Maple.

На практике испытательные частоты Oj = ¿%iax • (Oj f где щ вычисляются по формуле (12), могут задаваться с высокой, но ограниченной точностью. Поэтому для уменьшения вычислительной погрешности в разработанных процедурах Maple используются точные значения интерполяционных узлов, что улучшает работу алгоритма идентификации.

Выполнение равенств degA2N-i ~dsSA2N-i-\ =1 ' > = обеспечивает надежную работу алгорит-

мов КЧ и КЧА, так как это влечет определенность параметров {Ху, р^ }N 1 в формуле (9). Эти равенства в компьютерных вычислениях обычно выполняются автоматически, но возможность вырождения алгоритмов КЧ и КЧА позволило ввести в процедуры специальный численный параметр "tol" [10]. Удачный выбор пользователем параметра "tol" позволяет автоматически найти структурный параметр n при условии, что N > n . В разработанных процедурах при выборе tol = 0 алгоритмы заканчивают работу при

n = N .

Для идентификации исследуемого объекта в фиксированном диапазоне частот ,^max] вычисляется параметр у = / ^max и из априорных соображений выбирается число испытательных частот N.

Если выполнено неравенство у< COS((2N — Х)к / 4N) , то для идентификации следует использовать алгоритм КЧ, иначе, - алгоритм КЧА. По формуле (12), в которой для алгоритма КЧ у = 0 , находятся от-

,N ~

носительные частоты \G)jSfc—\ г применяемые для нахождения испытательных частот *

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

—1/2

Для сравнения трех методов рациональной аппроксимации действительной функции f (z) = z на

отрезке [у,1] , где у = 0,36 , в табл. 1 сведены результаты вычисления относительной погрешности 8 в % при использовании эллиптических функций (метод Золотарева) (4) - (6) с модулями k = 0,6, к' = 0,8,

алгоритма Паде-Чебышева и алгоритма Кронекера-Чебышева (N = n).

Таблица 1 Относительная погрешность 8 методов аппроксимации.

n Метод Золотарева Алгоритм Паде-Чебышева Алгоритм Кро- некера- Чебышева

2 0,651-10—2% 0,937 -10—2% 0,173 -10—1%

3 0,263 -10—4% —4 0,483-10 4% 0,114 -10—3%

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

—1/2

тельной области, успешно конкурируют с наилучшей рациональной аппроксимацией функции / (¿) = z ,

определяемой формулами (4) - (6). Наилучшая рациональная аппроксимация применяется, в частности,

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

Как показывают формулы (7) в задачах идентификации в частотной области необходима совместная аппроксимация АЧХ и ФЧХ.

—1/2 ~ г- “!

Для сравнения трех методов рациональной аппроксимации функции ИКр) = р , рЕУ]у^\, где

у — 0,36, к = 0,6, в табл. 2 сведены результаты вычисления относительной погрешности Зд в %, а в табл. 3 - Зф в % ^ = п) .

Таблица 2 Относительная погрешность Зд методов аппроксимации.

n Метод Золотарева Алгоритм Паде-Чебышева Алгоритм Кронекера-Чебышева Алгоритм Кронекера-Чебышева -Ахиезера

2 % 0 4 г- 7,53% % 8 3 2 -

3 % 4 4 1 1,48% - 0,411%

Таблица 3 Относительная погрешность 8ф методов аппроксимации.

n Метод Золотарева Алгоритм Паде-Чебышева Алгоритм Кронекера-Чебышева Алгоритм Кронекера-Чебышева -Ахиезера

2 9,27% 9,52% 0,572% -

3 % 1 7 1 1,71% - 0,313%

Лучшие результаты совместной аппроксимации АЧХ и ФЧХ с использованием алгоритмов КЧ и КЧА по сравнению с применением метода Золотарева и алгоритма ПЧ объясняются тем, что алгоритмы КЧ и КЧА

~ _______1/2

строят рациональную аппроксимацию функции Z(Jp) = p на множестве [—у,— ]у\ ^[уУ^]\ г состоящим из

двух отрезков мнимой оси.

—1/2

Рациональную аппроксимацию действительной функции / ^) = z на множестве 1-1-г№Л1 по-

строить не удается, так как продолжение этой функции на отрезок [—1, — у] в действительной области невозможно.

Следует также отметить, что при вычислении интерполяционных узлов {+/¿5^-}^ на множестве [—у, — ]у\ и [у/ау] параметр у непосредственно используется только в алгоритме КЧА в соответствии с

~ _______1/2

выражением (12). Это объясняет различие результатов аппроксимации функции г{р) = р 1,2 ' алгоритмами КЧ и КЧА, причем выбор одного из двух этих алгоритмов производится в зависимости от величины параметра у.

На рис. 1 и 2 приведены графики АЧХ и ФЧХ и их аппроксимации для иллюстрации особенностей соответственно методов Золотарева и КЧА при совместной аппроксимации АЧХ и ФЧХ с помощью рациональ-

~ ______1/2 ~

ной аппроксимации функции Z(^?) = р , /?е[уУ,у], где у = 0,36, к =0,6 , а N = п = Ъ .

Рис. 2. Аппроксимации АЧХ и ФЧХ с помощью алгоритма КЧА.

При N = п = 2 относительные частоты измерения (с пятью значащими цифрами) равны

Зх =0,38268, 32 =0,92388 , а при N = п = Ъ - 3^ = 0,43348 , 3^ =0,75153, 3^ =0,97041 . при N = п = 3 используется алгоритм КЧА, так как для у = 0,36 выполняется неравенство у> COS((2N — Х)л/(4^) = 0,25882

~ ~ ~ ~ 1/2 ~ ~ ~

Результат (р) рациональной аппроксимации модельной функции %(р) - Р г Р - ]3, 3 е [0,36; 1] с

помощью алгоритма КЧА, разложенный на простейшие дроби (используется пять значащих цифр), дается

соотношением

г3(р)-~

0,34430 0,40554 2,8003

- +-------------+ -

+ 0,083667 + 0,68253 ^ + 5,5359

Каждая простейшая дробь вида С 1/ (р + К 1С 1) определяет комплексное сопротивление КС-звена, образованного параллельным соединением резистора К и конденсатора С. Последовательное соединение различных КС-звеньев называется первой канонической формой Фостера, которая позволяет по извест-

ному разложению заданной функции импеданса на простейшие дроби построить ее реализацию в виде RC-двухполюсника [3].

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

Кроме известных канонических форм Фостера, Кауэра и Ли, автором получены еще две канонические формы. Также разработаны алгоритмы канонической реализации, позволяющие по заданной функции импеданса Zn (p) строить соответствующие канонические формы и различные их комбинации [11] . Эти алгоритмы реализованы в виде процедур в пакете Maple.

Каждый алгоритм канонической реализации устанавливает реализуемость функции Zn(Р) • Для этой цели самыми эффективными являются алгоритм Кауэра, соответствующий методу Стилтьеса проверки устойчивости заданного многочлена, а также модифицированный автором алгоритм Ли.

Так как при нормировке частоты 3=(0 1(0^^х справедливо равенство zni'P) = zn(P) - то получаем

рС = рС и С = Фщах^ • Поэтому реализация функции Zn(p), найденной в результате решения задачи аппроксимации, в виде RC-двухполюсника проблем не вызывает.

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

При нахождении RC-двухполюсника по результатам измерений значений АЧХ и ФЧХ следует учитывать погрешности измерений и дополнительные паразитные параметры, которые могут изменить вид RC-двухполюсника.

Для устранения этих эффектов, которые обусловлены наличием "шума", в алгоритмы КЧ и КЧА введен параметр "tol" [10]. Этот параметр подбирается пользователем из следующего соображения: требуется

найти такое минимальное n<N, при котором получаемая при аппроксимации функция Zn(p) , будет реализуемой. С использованием имитационного моделирования было проведено большое число экспериментов, подтверждающих результативность предложенного подхода для нахождения структурного параметра п.

7, Л -1/2

Так как в рассматриваемом модельном примере за счет применения известной функции Z(p) = p шум отсутствовал, то при подборе параметра "tol" получилось равенство tol = 0, которое повлекло соотношение N = n , означающее выбор структурного параметра п только из соображений точности. Полученное в результате имитационного моделирования равенство tol= 0 основано на возможности ап—1/2 —

проксимации функции f (z) = z аналитической функцией F ( z) класса S с любой степенью точности.

Таким образом, анализ модельного примера показывает преимущества алгоритмов Кронекера-Чебышева и Кронекера-Чебышева-Ахиезера для решения задачи структурно-параметрической идентификации в частотной области. Сравнение этих алгоритмов с методом конечно-частотной идентификации и взвешенным методом наименьших квадратов в частотной области рассмотрено в [10] и показано их применение к задачам автоматического регулирования.

ЛИТЕРАТУРА

1. Левин С. Ф. Теория измерительных задач идентификации // Измерительная техника. - 2001. -

№ 7. - С. 8-17.

2. Дамаскин Б. Б., Петрий О. А. Введение в электрохимическую кинетику. - М.: Высшая школа,

1975.

3. Балабанян Н. Синтез электрических цепей. - М.-Л.: Госэнергоиздат, 1961.

4. Крейн М. Г., Нудельман А. А. Проблема моментов Маркова и экстремальные задачи. - М.: Наука,

1973.

5. Бейтмен Г., Эрдейи А. Таблицы интегральных преобразований. Т. 2. - М.: Наука, 1970.

6. Ахиезер Н. И. Элементы теории эллиптических функций. - М.: Наука, 1970.

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

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

9. Бондаренко Л. Н. Метод Кронекера - Чебышева идентификации при ограниченных возмущениях // Труды Международной конференции "Проблемы автоматизации и управления в технических системах" -Пенза: Инф. изд. центр ПГУ, 2008. - С. 266-269.

10. Бондаренко Л. Н. Методы идентификации в частотной области при наличии шума // Известия ву-

зов. Поволжский регион. Технические науки. - 2009. - № 2. - С. 113-123.

11. Бондаренко Л. Н. Алгоритмические проблемы диагностики каскадных двухполюсников // Информационная измерительная техника. Труды университета. Межвуз. сб. науч. тр. - Вып. 32. - Пенза: Изд-во Пензенского государственного ун-та, 2008. - С. 3-12.