Конечно-частотная идентификация запаздывания с использованием фазовых сдвигов Текст научной статьи по специальности «Математика»

КОНЕЧНО-ЧАСТОТНАЯ ИДЕНТИФИКАЦИЯ ЗАПАЗДЫВАНИЯ С ИСПОЛЬЗОВАНИЕМ ФАЗОВЫХ СДВИГОВ1

Ливаткин П. А.2

(ООО "ТриаГруп", Москва)

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

Ключевые слова: теория автоматического управления, идентификация динамических объектов, идентификация запаздывания, конечно-частотная идентификация запаздывания.

1. Введение

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

1 Автор признателен д-ру физ.-мат. наук А.Г. Александрову за постановку задачи, а также коллективу лаборатории №7 Института проблем управления им. В.А. Трапезникова РАН, Москва, за активное обсуждение работы и полезные комментарии.

2Павел Анатольевич Ливаткин (pal2010@yandex.ru).

сываемых линейными дифференциальными уравнениями. Они различаются в зависимости от предположений о помехах измерения и внешних возмущениях, действующих на объект. Методы из первой группы рассматривают объекты, подверженные стохастическим воздействиям (например белый шум). Для их идентификации применяют методы стохастической аппроксимации и подходы, являющиеся различными модификациями метода наименьших квадратов [7]. Вторая группа методов рассматривает объекты с неизвестными ограниченными внешними возмущениями (при этом статистические характеристики полагают неизвестными). В этих случаях прибегают к конечно-частотной идентификации [9, 10] и рандомизированным алгоритмам [6].

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

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

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

длительности процессами идентификации на одной и той же положительной и отрицательной частоте.

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

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

В данной статье анализируется алгоритм конечно-частотной идентификации запаздывания из работы [1]. Далее будет показано, что он имеет ограниченные возможности для практического применения из-за особенностей используемого испытательного сигнала. Предлагается оригинальный алгоритм идентификации запаздывания на основе конечно-частотного подхода. Он основан на эквивалентности явления сдвига по фазе при прохождении гармоники через объект и её смещении в результате запаздывания. Это приводит к возможности модифицировать фильтры Фурье, введя отрицательный сдвиг по фазе в функцию свёртки. Результат расчета фильтров Фурье используется для восстановления коэффициентов передаточной функции объекта, которые сравниваются с результатом идентификации на другом наборе частот. Совпадение идентифицированных моделей говорит о том, что искусственно введённый сдвиг по фазе скомпенсировал запаздывание объекта, а значит равен ему. Рассматривается три разных способа оценки близости идентифицированных моделей.

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

2. Постановка задачи

Рассмотрим полностью управляемый, ассимптотически устойчивый объект с запаздыванием, описываемым уравнением

(1) (1^(1) +... + д,1№+ у® =

= кти(т\г - г) +... + кои(г - т) + / (г), где у(Ь) - выход объекта; и(Ь) - вход (управление); т - запаздывание, ограниченное сверху известным значением т* (т* > т > 0); кусочно-непрерывная функция /(Ь) - неизвестное ограниченное возмущение:

(2) |/ (*)| < /*,

где /* - некоторое неизвестное положительное число.

Необходимо решить задачу определения коэффициентов kj= 0, т; йг, г = 1,п, и запаздывания т объекта (1).

3. Анализ исходного алгоритма

В [1] решается задача определения коэффициентов kj, ] =0,т; йг, г = 1,п, и запаздывания т объекта (1).

Передаточная функция объекта (1) имеет вид

(3) ^ («) = № (8)е-тз,

т п

где W(з) = к(з)/й(з),к(з) = £ ^,й(з) = 1 + £

г=0 г=1

в - символ преобразования Лапласа при нулевых начальных условиях.

Для произвольной частоты 'Шг в [1] вводятся следующие обозначения частотных параметров:

рг = Ке(Ш и-шг)), & = 1т(Ш (¿-Шг)), () аг = Ее(Шг Ци,г)), & = 1т(Шг ^г)),

где ] - мнимая единица.

Для нахождения рг и фг потребуется вычислить коэффициенты Ш(в). Для этого в [1] применяется следующий подход. По определению:

(5) = аг +з!Зг,

¿(3 Ыг)

где ai и Pi вычисляются при помощи фильтров Фурье вида)

t^+ti tiF,+ti

(6) ai = J2- J y(t)sm(wit)dt,Pi = J y(t)cos(wit)dt,

где ¥ - длительность фильтрации, кратная целому числу периодов 2тг/ыг; £%р - момент начала фильтрации; у(Ъ) - выход объекта, к входному сигналу и(Ь) которого прибавили испытательный сигнал вида иег(Ь) согласно (7) (или сумму нескольких таких сигналов):

(7) ие&) = рг в'т (Ыг^.

Для объектов без запаздывания (Ш(,]ыг) = (.]ыг), так как т = 0. Значит и рг = аг,фг = @г) необходимо I = \(п + т + 1)/2\ уравнений конечно-частотной идентификации [1], позволяющих составить систему (8) и вычислить коэффициенты полиномов к(з) и ¿(в):

{Ш^г) = (рг + Зф{)с1(зы1), к(3^2) = (Р2 + ,

к(зыг) = (рI + зфг )й(зыг). Для объектов с запаздыванием ( = 0) система (8) примет

вид

( к(зи)г) = (р1 +зфх)е-Т^ги1й(зы{), (9) \ к(зм2) = (¥2 + 3Ф2)е-ТГШ2й(зы2),

\ к(зыг) = (р1 + зфг) е-ттй(зы1). К сожалению, неизвестен способ аналитически решить (9). В работе [1] для восстановления коэффициентов Ш(]ы) левая и

tL

F

правая части (9) домножаются на соответствующие комплексно-сопряженные выражения, вследствие чего пропадает сомножитель с запаздыванием, а уравнения принимают вид (10). Используя п + т + 1 частот, можно записать систему, решение которой даст коэффициенты полиномов k(s) и d(s) из алгебраических уравнений вида

(10) k(jwi)k(-jwi) e-Tjw. eTjw. = k{jwi)k(-jwi) =a¿ 2

d(jwi)d(-jwi) d(jwi)d(-jwi) г г'

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

= Re(W (X + jwx)),^x = Im(W (X + jwx)), ( ) ax = Re(Wz (X + jwx)),Px = Im(Wz (X + jwx)), где <£х,Фх вычисляются с помощью идентифицированной по формулам (6), (7) и (10) передаточной функции W (s); ax, fíx ищутся экспериментально при помощи испытательного сигнала вида uex(t) согласно (12) и (13):

(12) uex(t) = pxext sin (wxt),

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

(13), что будет пояснено ниже. При X > 0 длительность подачи идентифицирующего сигнала физически ограничена допустимой амплитудой входного сигнала, что также не позволяет получить асимптотическую сходимость фильтрации при действии внешнего возмущения.

В дальнейшем используются фильтры специального вида, являющиеся модификацией (6) для испытательного сигнала (12):

ax = ^ I е-Х*У (t)sin(wxt )dt,

+

(13) ^

Рх = 72х f e-Xty№ C°S (Wxt)dt,

где y(t) - выход объекта, к входному сигналу u(t) которого прибавили испытательный сигнал вида ue\(t) согласно (12) (или сумму нескольких таких сигналов); tx - длительность фильтрации

71

(целое число периодов sin (w*t)); tp - момент начала фильтрации.

Тогда верно следующее соотношение:

(14) k(X + jwx) е-r(x+3Wx) =

у ' d(X + jwx)

= (p* + jip\)e~rX(cos(wXT) - j sin(wXT)) = a* + jf3\. Из последнего равенства получим (14) получим т:

(15) r = ^ln ( Р* ' ( ) Т 2 Л m\a?x +Р1

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

Недостатки подхода:

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

2) подход можно использовать только для минимально-фазовых объектов;

3) отсутствие асимптотической сходимости процесса фильтрации (13) при наличии внешних возмущений.

Последний пункт вызван особенностью фильтрации при наличии внешнего возмущения. Рассмотрим фильтры (13) в условиях действия внешнего возмущения:

t$+tЛ

af\ = ^ J e~XtUf (t)sin(wxt )dt,

f Л

(16) 4+Л

Pfx = -Цр: / е xtyf (t)cos(wxt)dt,

где f( ) - составляющая выхода объекта, порожденная внешним возмущением.

При Л < 0 результат интегрирования растет быстрее чем , что приводит к всё более нарастающей ошибке. Данный

недостаток проиллюстрирован (см. рис. 1). На рисунке изображено изменение первого интеграла (16) при отсутствии испытательного сигнала при действии внешнего полигармонического возмущения (границы интегрирования 4 = о, ¿Л = 600). Явно видно, что ошибка нарастает. Минимум ошибки отмечен стрелкой.

5

3

1 о -1

О 100 200 300 400 500 В00

t, с

Рис. 1. График изменения фильтра Фурье af\ согласно формуле (16)

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

4. Конечно-частотная идентификация объекта с запаздыванием с использованием фазовых сдвигов

4.1. МАТЕМАТИЧЕСКОЕ ОБОСНОВАНИЕ

Вернёмся к уравнениям (6). Запись в них измеряемого сигнала как y(t) осложняет анализ запаздывания. Поэтому представим его как (17) для объекта без запаздывания и как (18) для объекта с запаздыванием в г:

(17) y(t) = A(jwi) sin (wit + S(jWi)),

(18) y(t) = A(jwi) sin (wit - WiT + S(jWi)),

где A(jwi) = pi ■ \ W(jwi)\ - изменение амплитуды pi входного сигнала u(t) начастоте w, объектом W(s); ó(jw,¿) = arg[W(jwj,)] -изменение объектом W(s) фазы входного сигнала u(t) начастоте wi

Хотелось бы ввести в уравнение (18) параметр в так, чтобы выполнились следующие требования:

1) компенсировался вводимый запаздыванием член w,t, приводя уравнение (18) к виду (17), для которого уже разработана теория конечно-частотной идентификации;

2) параметр имел математическую связь с запаздыванием, позволяя по в вычислить т;

3) единственное значение параметра в обеспечивало компенсацию WiT.

С учётом требований уравнения (6) дополняются (18) и искусственно введённым сдвигом по фазе Wi6:

t%+ti

ai,e = f A(jwi)sin(wit-wiT+wi9+S(jwi))sin(wit)dt, (19) ' "

I%,e = jr^i f A(jwi)sm(wit-WiT+Wi6+S(jwi)) cos (wit)dt.

* ^

Уравнения (19) можно получить на практике несколькими путями. Например, вводя w,9 в испытательный сигнал (7), вычисляя ai:e и Pj,,^, сдвигая в до тех пор пока запаздывание не будет скомпенсировано:

(20) Uei(t) = pi Sin (Wit + Wiff),

y(t) = A(jwi) sin (wit - WiT + Wiв + 5(jwi)).

Уравнения (19) и (20) хорошо показывают, на чём основан предлагаемый метод идентификации объекта с запаздыванием, однако при использовании на практике такой подход будет занимать много времени, либо памяти, либо вычислительных ресурсов. Поэтому в тексте работы будут встречаться отсылки к (19), поясняющие алгоритм, хотя применимы на практике эквивалентные им уравнения (34), нуждающиеся в единственном экспериментальном вычислении параметров ai и Д.

При в = т результат вычисления фильтров Фурье (19) становится эквивалентен получению фильтров Фурье от объекта без

запаздывания W(s) (согласно (4)) на той же частоте:

éF+ti

Pi = f A(jwi) sin (wit + S(jwi)) sin (wit)dt

(21)

Pit1

tFp+f

2

fi

F

= P2fi f A(jw¿) sin (wit + S(jwi)) cos (wit)dt.

i

F

Из уравнения (19) видно, что при изменении введённого фазового сдвига в величины (ai,e, fi,e) будут повторяться через период Pi = W. Следовательно, не только aiT = pi и fiT = éi, но

Wí ' '

и

(22) aifi = pi,6 = т + P-n, n eN,

( ) Pi,e = Ф-,д = т + Pin, n eN.

Частотные параметры ( ai,e ,fite) используются для записи системы уравнений конечно-частотной идентификации. Решение системы (23), аналогичной (8), из I = \(n + m + 1)/2\ линейных уравнений относительно коэффициентов полиномов k(s) и d(s) соответствует фиктивной передаточной функции Wg (s), для которой, очевидно, выполняется соотношение ai,e +jfi,e = We (jwi):

Ík(jwi) = (ai,e + jfi,e )d(jw\), k(jw2) = (a2,o + jf2,o )d(jw2),

k(jwi) = (at,e + jfi,e )d(jw¡). Каждое из уравнений данной системы (23) будет обладать свойством (22), но при этом иметь свой уникальный период P-.

Введём понятие наименьшего общего кратного (Icm) для рациональных чисел а,Ь,с,...: 1 cm(a, b,c,...) - это наименьшее рациональное число, которое при делении на любой из аргументов даёт целое число. Наименьшее общее кратное с иррациональным числом даёт те.

Тогда решение системы (23), являющееся коэффициентами передаточной функции We (s), будет совпадать с решением (8) при

(24) 0 = T + P (We),

где P(W$) - период, с которым повторяется величина всех коэффициентов (ai,e, fi,e) системы уравнений (23) согласно (25):

(25) _ Р (Шв ) = 1ст(Р1,Р2,...,Р1),

где Рг,г = 1,1, - периоды на частотах из системы (23).

Решение системы (8) неизвестно. Чтобы его отыскать, необходимо подать дополнительный дополнительные испытательные сигналы (хотя бы 1, что является преимуществом по сравнению с методом, предлагаемым в [1]) на частотах (г = I + 1,1 + р). Точное количество выбирается исходя из неравенства (26). Таким образом могут быть получены ещё р ^ 1 уравнений для построения систем вида (23) из I уравнений. Благодаря этому на разных наборах частот можно составить до М = вариантов систем. Для обозначения передаточных функций, полученных решением таких систем, введём обозначение (в),д = 1,М. Поскольку

(5) = (в) = ... = Wq,т(в) = Ш(в), то равенство между собой идентифицированных (в),д = 1,М, может служить критерием обнаружения истинного запаздывания: в = т.

Выберем частоты ,г = 1,1 + р, так, чтобы период повтора всех систем частотных уравнений составил

(26) Рт = 1ст(Р ),Р (Ш2,в),..., Р (Шм,в)) > г *.

Из этого следует, что идентифицированные на разных наборах частот передаточные функции (в),д = 1,М, совпадут между собой единственный раз в допустимом диапазоне. Это условие позволит однозначно идентифицировать запаздывание, и, как следствие, рассчитать коэффициенты передаточной функции объекта без запаздывания:

(27) 0 < в = т<Т*.

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

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

4.2. ВАРИАНТЫ ОЦЕНКИ БЛИЗОСТИ ПЕРЕДАТОЧНЫХ ФУНКЦИЙ

1. Близость корней полиномов.

Обозначим корни полиномов числителей kq,g (s) передаточных функций Wq,e (s) как Rq,e = [г\,q,e, r2,q,0,..., rm,qß ], а корни знаменателей dq>e (s ) - как Xq,d = [xi>q>e, x2>q>e,..., xn>q>e ]. При в = т будут выполняться равенства

(28) Ri,e = Ri+\,e; хг,в = хг+\,в; г = l,M — 1. Поскольку при экспериментальном определении (щ, ßi)

присутствует некоторая ошибка, то следует искать в, минимизирующую функционал М-1 м

(29) Ег = min ^ ^ (R3,e - R%e)(Rjfi - Ri,e)H+

j=l i=j+\

+ (Xj,e — Xi,e )(Xj,e — Xi,e)H. Минимизация этого функционала требует «правильной» упорядоченности корней в векторах Riß и Xiß, что теоретически является вычислительно трудной задачей. Однако если корни полиномов передаточных функций близки друг к другу либо порядок полиномов объекта не превышает 2, то правильное упорядочивание является простым.

2. Близость коэффициентов полиномов. Альтернативной, более простой с точки зрения вычислений,

но менее точной по результатам, может служить минимизация суммарного «расстояния» между коэффициентами идентифицированных передаточных функций Wqß (s). При равенстве коэффициентов полиномов их корни также будут равны, что приведёт к выполнению (28), а значит, эквивалентно (29). При этом (30) не требует вычисления корней полиномов и правильной их упорядоченности, что представляло сложности в первом способе. М-1 м

(30) Е2 = min V V (К3,в - Ki,g)(Kj,g - Ki,g)т+

j=l i=j+\

+ (Djt0 - Di,0)(Dj,0 - Di,0)т,

где К^д,Кг,в - векторы коэффициенты полиномов к^д(я) и кг,в(я) соответственно; Dj,вв - векторы коэффициенты полиномов ($) и йгв (в) соответственно.

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

3. Близость измеренных и расчётных частотных параметров. Третьим вариантом является подстановка некоторой новой частоты в передаточную функцию (^г), которая была получена без использования -шх. То есть сравнение расчётного значения Ш,,ваых) с измеренными значениями (19) по расстоянию на комплексной плоскости:

(31) 1<Ух,в+з@х,в - а^х) =

= Цых)} - ах,в)2 + (1т{Ш,,в^х)} - рх,в)2.

Выражение (31) должно равняться 0 при в = т, так как точка (ах,т,Рх,т) принадлежит годографу Шт(в).

Поскольку при экспериментальном определении (ах,в,@х,в) присутствует некоторая ошибка, то следует искать в, минимизирующую функционал

м р

(32) Ез = ппп^ ^ 1а.х,в + зРх,в - (jWх)|.

,= 1х=1

4.3. УПРОЩЕНИЕ ВЫЧИСЛИТЕЛЬНОЙ СОСТАВЛЯЮЩЕЙ МИНИМИЗАЦИИ ФУНКЦИОНАЛОВ Во всех трех случаях необходимо проверить большое количество значений в, что приводит к повторному вычислению интегралов (19). Для ускорения поиска в получим формулы для вычисления а^ в и @г,в с меньшей вычислительной сложностью. Идея заключается в том, что экспериментально получаемые параметры аг и Д идентифицируемого объекта могут быть преоб-

разованы к параметрам аналогичного объекта Шв (,]шг) с произвольным сдвигом по фазе в подобно введению запаздывания:

(33) Шв(зыг) = Ш2(зыг)ев™ = (он + з!Зг)ев™ = аг,в + зРг,в.

Применяя формулу Эйлера к ев:>ш* в (33), получим следующие равенства, позволяющие быстро вычислить ( аг,в, Рг,в) по экспериментально получаемым по формулам (6) значениям (0Ог,Рг):

(34) Огв = Ог СОв(Шгв) - Р вШ^в), Рг,в = Рг в) + Ог в).

4.4. ОБЩИЕ РЕКОМЕНДАЦИИ ПО ВЫБОРУ ЧАСТОТ

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

1. Не стоит выбирать их слишком близкими. На близких частотах будут получаться близкие значения (Ог,Рг), что может приводить к некорректным решениям системы (23) из-за близости к вырожденности самой системы.

2. Иррациональные числа кажутся хорошим выбором, поскольку могут обеспечить РТ = те при добавлении всего одной дополнительной испытательной частоты (р = 1). Благодаря этому отпадает необходимость заранее знать верхнюю оценку запаздывания. К сожалению, их невозможно реализовать физически, тем не менее использование их приближений может существенно увеличить Рт.

Алгоритм 1.

1. Зададимся в = 0, т* и шагом А в (чем меньше, тем лучше, исходя из вычислительных мощностей).

2. Зададимся набором частот шг,г = 1,1 + р, удовлетворяющим критерию (26). Для каждой из них экспериментально найдем значения ( аг, Рг) согласно (6).

3. По формуле (34) с учётом текущего значения в вычислим частотные параметры (ai,e,pi,e).

4. Выберем до М = ^^ур наборов значений (ai,$, pi,e), исключая р пар фазочастотных параметров. Для каждого набора составим систему из I уравнений вида (23). Решая её найдем коэффициенты Wi,8 (s),i = 1,М.

5. Оценим близость Е полученных передаточных функций Wi,e(s),i = 1,М, по одному из предложенных критериев (29), (30), (32) или их сумме.

6. Если Е = 0, т.е. полученные Wi,e (s), i = 1, М, совпадают, то в = т согласно (24) и задача идентификации объекта с запаздыванием решена. Конец алгоритма 1.

7. Если в < т*, то сохраним текущую пару (в, Е). Будем увеличивать в на величину АО и возвращаться к пункту (3) после каждого изменения в.

8. Среди сохраненных пар (в, Е) ближайшим к т значением в будет то, у которого значение критерия Е минимально. Конец алгоритма 1.

5. Примеры

5.1. РЕЗУЛЬТАТЫ ИСХОДНОГО АЛГОРИТМА

В [1] приведён объект с запаздыванием и ограниченным возмущением

(35) 0,7ij(t) + 0,8y(t) + y(t) = 0,4ñ(t - 3) + u(t - 3) + f (t), где внешнее возмущение

(36) f (t) =2sign[sin(5t)].

В тексте статьи предлагалась двухэтапная идентификация, где на первом этапе к входу объекта прибавляется испытательный сигнал вида (7) длительностью около 186 с (21 период наименьшей испытательной частоты) с интервалом дискретизации по времени h = 0,01:

(37) uei(t) = 0,05 sin (0,707t) + 0,075 sin (1,41t) +

+ 0,1 sin (2,12t) + 0,15 sin (2,83t).

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

(38) ueX(t) = 0,1 e-0'01t sin (0,707t).

В результате идентификации был получен объект:

(39) 0,65y(t) + 0,83y(t) + y(t) = 0,32U(t - 3,48) +

+ 1,06u,(t - 3,48) +f(t).

5.2. РЕЗУЛЬТАТЫ ФАЗО-ЧАСТОТНОИ ИДЕНТИФИКАЦИИ

С АНАЛОГИЧНЫМ ИСПЫТАТЕЛЬНЫМ СИГНАЛОМ

Для сравнения эффективности алгоритмов ко входу объекта был прибавлен испытательный сигнал близкий к (37): с теми же частотами и амплитудами компонент, но меньшим числом гармоник. Согласно алгоритму 1:

1. Выбраны параметры т* = 10 и А в = 0,001.

2. Выбраны две частоты w1 = 0,707 и w2 = 1,41 (I = 2), а также одна дополнительная частота w3 = 2,12 (р = 1). На них проведена идентификация длительностью в 21 период наименьшей испытательной частоты. Интервал дискретизации по времени h = 0,01. Испытательный сигнал принял вид

(40) uei(t) = 0,05 sin (0,707t) + 0,075 sin (1,41t) +

+ 0,1 sin (2,12t).

Получены пары значений: (a^ f1) = (-0,9723; -0,7862), (a2]p2) = (0,7988,0,6021) и (a3; fi) = (-0,1392; -0,452).

3. По формуле (34) с учётом текущего значения в они пересчитаны в частотные параметры ( a1to;f1te), (a2,$;f2,e) и (a3,e ;Рз,в).

4. Из полученных частотных параметров создаётся М = 3 системы уравнений вида (23):

k(jw1) = (a1,e + 3 fifi )d(jwi),

k(jw2) = (a2,e + jfh,e )d(jw2);

k( w1) = (a1,e + 3f1.fi )d(jwi),

k(jw3) = (аз,в + 3 fifi )d(jw3);

k( w2) = (a2,e + jf2fi )d(jw2),

k(jw3) = (аз,в + 3 fsfi )d(jw3).

5. Вычислены значения Е1, Е2 и Е% полученных передаточных функций Ш1>в(з),Ш2,в(з),Ш3,>в(в) по критериям (29), (30), (32) соответственно.

6. Фаза сдвигалась на величину АО: в = в + АО и повторялись пункты с (3) по (5) до тех пор, пока 9 < т* = 10. Получаемые на каждом цикле Е сохранялись в соответствующие массивы. Полученные зависимости Е(в) приведены в соответствии с графиком на рис. 2.

01 23456789 10 6

01 23456783 10

е

01 23456789 10

е

Рис. 2. График зависимости функционалов (29), (30) и (32) от в при испытательном сигнале из [1] 7. Выбраны значения в, соответствующие минимумам Е(в). Результаты идентификации с использованием каждого из предложенных критериев (29), (30) или (32):

а) по критерию близости корней полиномов (^1):

(42) 0,788у(г) + 0,739у(Ь) + у(Ь) = 0,5455и(г - 3,054)+

+ 0,8827и(Ь - 3,054) + / (г).

б) по критерию близости коэффициентов полиномов (^2):

(43) 0,7184у^) + 0,7036у(Ь) + у(Ь) = 0,3981й(г - 2,976) +

+ 0,9509и(г - 2,976) + /(Ь).

в) по критерию близости измеренных и расчётных частотных параметров (Ез):

(44) 0,7184у^) + 0,7036у(Ь) + у(Ь) = 0,3981й(г - 2,976) +

+ 0,9509и(г - 2,976) + /(Ь).

5.3. РЕЗУЛЬТАТЫ ФАЗО-ЧАСТОТНОЙ ИДЕНТИФИКАЦИИ С ПОДОБРАННЫМИ ЧАСТОТАМИ

Для демонстрации эффективности алгоритма произведено моделирование процесса идентификации объекта из [1] согласно алгоритму (1) с использованием каждого из предложенных критериев:

1. Выбраны параметры т* = 10 и АО = 0,001.

2. Исходя из верхней оценки запаздывания т* и (26), учитывая [11], были выбраны две частоты w1 = 0,2ж и = 0,8^ (I = 2), а также одна дополнительная частота w3 = ж (р = 1). На них проведена идентификация длительностью в 80 с. Интервал дискретизации по времени к = 0,01.

(45) иег(г) = 0,05ё\п(0,2^г)+ 0,08 ё\п(0,8^г)+ 0,1 ё\п(^г). На выбранных частотах получены пары значений: (а.1; Д) = = (-0,4828;-1,1647), (а.2; Р2) = (-0,3747; -0,0207) и (аз; Рз) = (0,0662; 0,2493).

3. По формуле (34) с учётом текущего значения 0 они пересчитаны в частотные параметры (а1,в; Р1,в), (а2,в; Р2,в) и (аз,в; Рз,в).

4. Из полученных частотных параметров создаётся М = 3 системы уравнений, аналогично (41).

5. Вычислены значения Е1, Е2 и Ез полученных передаточных функций Ш1>в(з),Ш2,в(8),Шз,>в(в) по критериям (29), (30), (32) соответственно.

6. Фаза сдвигалась на величину АО: 0 = 0 + АО и повторялись пункты с (3) по (5) до тех пор, пока 0 < т*. Получаемые

на каждом цикле Е(в) сохранялись в соответствующие массивы и приведены в соответствии с графиком на рис. 3.

7. Выбраны значения в, соответствующие минимумам Е(в). Результаты идентификации с использованием каждого из предложенных критериев (29), (30) или (32):

а) по критерию близости корней полиномов (Е1):

(46) 0,8303у(Ь) + 0,5313у(Ь) + у(Ь) = 0,5516й(г - 3,006) +

+ 0,8694и,(Ь - 3,006) + /(г).

0 1 23456739 10

6

01 23456733 10

0

01 23456739 10

е

Рис. 3. График зависимости функционалов (29), (30) и (32) от в

б) по критерию близости коэффициентов полиномов (Е2): (47) 0,784у(Ь) + 0,505у(Ь) + у(Ь) = 0,4985и(г - 2,987)+

+ 0,8966и,(Ь - 2,987) + /(г).

в) по критерию близости измеренных и расчётных частот-

ных параметров (Ез):

(48) 0,7817у(г) + 0,5037у(Ь) + у(Ь) = 0,4958й(г - 2,986) +

+ 0,8979и(Ь - 2,986) + /(Ь).

6. Результаты

Можно видеть существенное улучшение точности идентификации запаздывания. Погрешность в определении т, рассчитанная по формуле 100(1 - в/т), снизилась с 16% в (39) до 0,2% для первого критерия и до 0,43% и 0,47% - для второго и третьего соответственно. Эти результаты получены при вдвое меньшей длительности идентификации и сократившемся количестве испытательных гармоник!

Для оценки поведения различных критериев оценки близости приведен график (см. рис. 3), полученный в ходе фазо-частотной идентификации с оптимальными частотами. На нём видно, как все функционалы стремятся к нулю в окрестности т. При этом на некоторых кривых имеются локальные минимумы, которые при негативном стечении факторов (недостаточно долгая идентификация, внешнее возмущение большой мощности, близость испытательной частоты и частоты внешнего возмущения) могут быть приняты за истинное запаздывание.

Анализируя полученные кривые, можно прийти к выводу, что оценка с использованием первого критерия (чётко выраженный минимум) менее подвержена влиянию внешних возмущений. На графиках функционалов (30) и (32) видны локальные минимумы, близкие к глобальному. Сумма критериев на одинаковых в, по всей видимости, позволит отбросить большинство таких локальных минимумов. Погрешность в определении т ниже у первого способа.

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

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

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

7. Выводы и перспективы

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

Приведённый в работе алгоритм был успешно применен на практике в системе автоматизированного управления вакуумной дуговой печью на АО «Металлургический завод "Электросталь"».

Литература

1. АЛЕКСАНДРОВ А.Г., ОРЛОВ Ю.Ф., ПАЛЕНОВ М.В.

Конечно-частотная идентификация объектов с запаздыванием // Автоматика и телемеханика. - 2014. - №2. - С. 5-15.

2. АЛЕКСАНДРОВ А.Г., ОРЛОВ Ю.Ф. Сравнение двух методов идентификации при неизвестных ограниченных возмущениях // Автоматика и телемеханика. - 2005. - Т. 66, №10. -С. 128-147.

3. АЛЕКСАНДРОВ А.Г. Частотные регуляторы // Автоматика и телемеханика. - 1991. - №1. - С. 22-33.

4. АЛЕКСАНДРОВ А.Г. Конечно-частотная идентификация: границы частот испытательного сигнала // Автоматика и телемеханика. - 2001. - №11. - С. 3-14.

5. АЛЕКСАНДРОВ А.Г. Метод частотных параметров // Автоматика и телемеханика. - 1991. - №12. - С. 3-15.

6. ГРАНИЧИН О.Н., ПОЛЯК Б.Т. Рандомизированные алгоритмы оценивания и оптимизации при почти произвольных помехах. - М.: Наука, 2003.

7. ЛЬЮНГ Л. Идентификация систем. - М.: Наука, 1991.

8. AHMED S., HUANG B., SHAH S. Process identification from sinusoidal test data by estimating step responce // Preprints of the 15th IFAC Symposium on System Identification. - Saint-Malo, France, 2009. - P. 396-401.

9. ALEXANDROV A.G. Finite-frequency method of identification // Proc. of 10th IFAC Symposium on System Identification. Preprints. - Vol.2 - Copengagen, Denmark, 1994. - P. 523-527.

10. ALEXANDROV A.G. Finite-frequency identification and model validation of stable plant // Proc. of 14th World Congress of IFAC. Preprints. - Vol. H - Beijing, China, 1999. -P. 295-301.

11. ALEXANDROV A.G. Finite-frequency identification and model validation of stable plant // Proc. of 16th World Congress of IFAC. Preprints. - Prague, Czech Republic, 2005.

12. ALEXANDROV A., ORLOV YU., PALENOV M. Finite-frequency identification of plant with time delay // Preprints of the 16th IFAC Simposium on System Identification. - Brussels, Belgium, 2012. - P. 66-70.

13. BJORKLUND S., LJUNG L. A review of time-delay estimation techniques // Proc. of the 42nd IEEE Conference on Decision and Control. - Hawaii, USA, 2013.

14. M'SAAD M., FARZA M. Identification of continuous-time linear systems with time-delay // Preprints of the 15th IFAC Symposium on System Identification. - Saint-Malo, France, 2009. - P. 898-903.

15. WONG K.Y., POLAK E. Identification of linear discrete time systems using instrumental variable method // IEEE Trans. Automat. Control. - 1967. - Vol. AC-12. - - P. 707-718.

16. YAMAGUCHI H., NAGUMO K., ABE N. Iterative design using IMC for input time-delay systems with disturbance observer // Preprints of the 18th IFAC World Congress. -Milano, Italy, 2011. - P. 2430-2435.

FINITE-FREQUENCY IDENTIFICATION OF A DELAY VIA SLIDING PHASE LAG

Pavel Livatkin, TRIA GROUP, Moscow (pal2010@yandex.ru).

Abstract: This article proposes a phase sliding improvement of finite-frequency identification algorithm for a linear stable plant with time-delay in presence of unknown-but-bounded external disturbances (with unknown stochastic characteristics). Finite-frequency identification algorithm feeds to the plant's input a testing signal, that consist of a single harmonic or a sum of them. Phase sliding improvement allows to identify two unknown values, such as coefficients of the plant's transfer function or time-delay value for the successful identification and more than one harmonic if it is needed to increase accuracy or identify a long time-delay value. There are two ideas in the article. Firstly, phase sliding for each harmonic may be opposite for phase deviation, caused by plant's time-delay. Time-delay, phase sliding and harmonic's frequency are analytically related. Secondly, plant's transfer function will be the same for different combinations ofharmonics ifthey are opposite to plant's time-delay. For a definite solution of the delay identification problem the upper bound of the time-delay value is required. The article describes how to select theoretically optimal harmonics, providing time-delay identification without known upper bound. There are three algorithms for the detection equality of phase-slided and time-delayed values, which differ one from each other computational complexity and sensitivity to external disturbances.

Keywords: control systems engineering, dynamical systems identification, delay identification, finite-frequency delay identification.

УДК 681.5 ББК 32.966

DOI: https://doi.org/10.25728/ubs.2020.84.4

Статья представлена к публикации членом редакционной коллегии П.С. Щербаковым.

Поступила в редакцию 15.10.2019. Дата опубликования 31.03.2020.