Лекции 2014 и 2010 годы (Лекции по курсу Системы измерения и оценки параметров)

Посмотреть архив целиком

Глава 4. Идентификация и измерения

Боевкин В.И.

БОЕВКИН ВИКТОР ИВАНОВИЧ Лекции по курсу СИСТЕМЫ ИЗМЕРЕНИЯ И ОЦЕНКИ ПАРАМЕТРОВ И ХАРАКТЕРИСТИК. 9 семестр




1. Введение

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

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

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

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

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



2. Разложение процесса по произвольному (неортогональному) базису

2.1. Алгоритм неортогонального разложения

При обработке экспериментальных функций времени может возникнуть необходимость представления этих функций в виде аналитических выражений. Непрерывный процесс f(t), определенный на интервале T, можно разложить по заранее выбранной системе базисных функций ф1(t)…фn(t), где n – число базисных функций. Разложение S(t) имеет вид:

S(t) = Sk=1n(Ck фk ) ; f(t) = S(t) + d(t) (2.1)

Здесь

d(t) = f(t) – S(t) - ошибка разложения (2.2)

Введем понятие энергии Ey функции (сигнала) y(t):

Ey = [0Ty2(t)]dt (2.3)

В качестве критерия, характеризующего удаление S(t) от f(t) на интервале T, то-есть характеризующего качество разложения, примем энергию ошибки:

Ed = [0Td2(t)]dt (2.4)

При заданных f(t) , T и фi(t) энергия ошибки является функцией от коэффициентов Ci разложения. Коэффициенты разложения будем выбирать из условия минимума энергии ошибки, для чего приравняем нулю частные производные

DEd/DCi = 0 ; i = 1…n . (2.5)

Из (4) и (2) можно получить:

Ed = [0T{f(t) – S(t)}2dt = Ef – A + Es ; (2.6)

Ef = [0T f 2 (t)] dt , A = -2 [0T f (t) S(t)] dt , Es = [0T S 2 (t)] dt

С учетом (6) уравнения (5) станут:

DEd/DCi = DEf/DCi - DA/DCi + DEs/DCi = 0; i = 1…n (2.7)

Используя (1) и (6) , выпишем частные производные:

DEf/DCi = 0; DA/DCi = 2[0T f(t) DS/DCi dt (2.8)

DEs/DCi = 2 [0T S(t) DS/DCi ]; DS/DCi = фi (t)

Подставив (8) в (7), получим:

Sk=1n (Ck Ui,k ) = Vi; i = 1…n . (2.9)

Здесь введены обозначения:

Ui,k = [0T фi (t) фk (t)] dt; Vi = [0T фi (t) f (t)] dt. (2.10)



Линейная алгебраическая неоднородная система уравнений (9) вместе с обозначениями (10) определяет значения коэффициентов разложения Ci , при которых энергия Ed ошибки разложения минимальна. Используя обозначения (10), энергию ошибки (4) при произвольных значениях коэффициентов Ci можно представить:

Ed = Ef – 2 Si=1n(Ci Vi) + Si=1n{Ci Sk=1n (Ck Ui,k)} . (2.11)

При оптимальных значениях коэффициентов, вычисленных из (9), энергия ошибки (11) имеет вид:

Ed = Ef - Si=1n(Ci Vi) = Ef –Es . (2.12)

Приведенные выше результаты имеют достаточно прозрачное геометрическое представление. Введем в рассмотрение пространство H (бесконечномерное), элементами которого являются функции, определенные на интервале T, и удовлетворяющие условиям Дирихле, то – есть ограниченные и имеющие конечное число экстремумов и разрывов первого рода. Любая такая функция представляется в пространстве H точкой или, что эквивалентно, вектором. Определим скалярное произведение двух векторов в пространстве H следующим образом:

^x(t) , y(t)^ = [0T x(t) y(t)] dt (2.13)

Пространства, в которых определена операция скалярного произведения, называются Гильбертовыми. Они обладают геометрическими аналогиями, вследствие которых суждения, полученные из геометрического рассмотрения, имеют силу доказательств. Квадрат модуля любого сигнала y(t) в пространстве H равен энергии этого сигнала, что видно из определений энергии (3) и cкалярного произведения (13):

^y(t) , y(t)^ = [y(t)]2 = [0T y2(t) dt = Ey (2.14)

Векторы y(t) и x(t) ортогональны, если их скалярное произведение равно нулю:

^x(t) , y(t)^ = [0T x(t) y(t)] dt = 0 (2.15)

Базисные функции фi(t) изображаются в пространстве H векторами, которые формируют гиперплоскость Sn размерностью n, то есть набор функций фi (t) является неортогональным базисом гиперплоскости Sn. Разложение S из (1) всегда лежит в гиперплоскости Sn а коэффициенты Ci являются координатами S в неортогональном базисе гиперплоскости Sn. Разлагаемая функция f(t), в общем случае, не лежит в гиперплоскости Sn. Ошибка разложения d(t) изображается в пространстве H вектором, соединяющим концы векторов f(t) и S(t). Квадрат длины вектора ошибки, согласно (14), равен энергии ошибки, т.е. принятому критерию близости (4). Из геометрических соображений ясно, что длина [d] вектора ошибки (а вместе с ней и энергия ошибки) будет минимальной, если вектор ошибки d(t) ортогонален гиперплоскости Sn. Последнее выполняется, когда вектор ошибки d(t) ортогонален ко всем базисным функциям фi(t), формирующим гиперплоскость Sn. Принимая во внимание определение скалярного произведения (13) и условие ортогональности (15), можно записать:

^d(t) , фi(t)^ = [0T d(t) фi(t)]dt = 0

Легко видеть, что подставив сюда d(t) = f(t) – S(t) из (2) и используя обозначения (10), получим систему уравнений (9). Величины (10) при этом являются скалярными произведениями:

Ui,k = ^фi(t) , фk(t)^ = [0T фi (t) фk (t) dt ; (2.16)

Vi = ^f(t) , фi(t)^ = [0T f(t) фi (t) dt .

Матричная запись системы уравнений (9) имеет вид:

U*C = V (2.17)

Здесь U – квадратная матрица размерностью n*n с элементами Ui,k , V – вектор правых частей размерностью n*1 с элементами Vi , C– вектор неизвестных размерностью n*1 с элементами Ck..

4.2.2. О сходимости разложений



При заданном базисе фi(t) коэффициенты Ui,k системы уравнений (9) могут быть вычислены заранее. В этом случае определение коэффициентов разложения Ci сводится к вычислению правых частей Vi по (16) и решению системы линейных уравнений (9) или, что эквивалентно, матричного уравнения (17). При большой размерности n базиса последняя операция достаточно трудоемка. Если выбрать такие базисные функции, которые в пространстве H ортогональны между собой попарно, то есть

i(t) , фj(t)^ = [0T фi(t) фj(t)]dt = 0; (2.18)

i = 1…n; j = 1…n ; i [=] j ,

матрица U становится диагональной, а система (9) – разрешенной. Разложение по ортогональному базису при возможности увеличения размерности базиса до бесконечности и выполнении некоторых других условий называется обобщенным рядом Фурье [Л2]. Простота решения системы уравнений (9) в большинстве случаев определяет использование ортогональных базисов (обобщенных рядов Фурье), в том числе и тригонометрических. Однако, в некоторых задачах предпочтение может быть отдано неортогональным разложениям. Степень гладкости разлагаемой функции можно характеризовать порядком R младшей производной, в которой на интервале разложения существуют разрывы первого рода (конечные разрывы). Известно [Л3], что сходимость тригонометрического ряда Фурье, под которой понимается уменьшение коэффициентов ряда с ростом его номера, имеет порядок 1/nR. Однако, при сколь угодно высокой степени гладкости разлагаемой функции внутри интервала разложения, несовпадение значений разлагаемой функции в начале и конце интервала разложения эквивалентно разрыву в функции, что приводит к весьма слабой сходимости ряда порядка 1/n. Это является следствием периодичности ортогонального тригонометрического базиса на интервале разложения. Возможность улучшения сходимости за счет отказа от ортогональности базиса проиллюстрируем простым примером. Разложим функцию f(t) = exp(-t) на интервале (0…2) по базису:



1, cos(pxt), sin(pxt), cos(2pxt),sin(2pxt), … cos(qpxt),sin(qpxt) (4.2.19)

Здесь q – число гармоник в базисе. Базис (19) ортогонален при x = 1. Отличие величины x от единицы характеризует удаление базиса от условий ортогональности. На рис.4.1 представлена зависимость от x среднеквадратической ошибки разложения

для базисов различной размерности: q = 1, q = 2, q = 3.



На Рис.4.2 представлены графики разлагаемой функции разложения по ортогональным базисам различных размерностей (кривые q = 1, q =2, q = 3) и по неортогональному базису размерности q = 1 (кривая x = 0.1). Из графиков, представленных на рисунках 4.1 и 4.2, можно видеть, что увеличение размерности ортогонального базиса не приводит к существенному уменьшению ошибки разложения. Применение же неортогонального базиса позволяет существенно уменьшить ошибку разложения уже при малой размерности базиса. Увеличение объема расчетов на решение уравнений (17) может компенсироваться увеличением точности разложения.

4.2.3. Разложение дискретных функций

Экспериментальные данные обычно являются дискретными функциями времени, то есть представляют собой набор отсчетов:

f(tj) = f(j); j = 1…N.



Для разложения дискретных функций базисные функции и само разложение также должны быть дискретными:

фi(tj) = фi(j), S(j) = Si=1 n [Ci фi(j)], i = 1…n, j = 1…N

Функциональное пространство H, введенное выше, становится N – мерным, а скалярное произведение определяется соотношением:

^x(j), y(j)^ = Sj=1N x(j) y(j)

Энергия ошибки (4) (минимизируемый критерий) имеет вид:

Ed = Sj=1Nd2(j) = Sj=1N(f(j) – S(j))2

Уравнения (9) и (17) для определения оптимальных коэффициентов разложения остаются в силе, но в параметрах Ui,k и Vi интегралы заменяются суммами:

Ui,k = Sj=1Ni(j) фk(j)]; Vi = Sj=1Ni(j) f(j)] (4.2.20)



4.2.4. Фильтрующие свойства разложений

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

f(j) = S0(j) + ah(j); S0(j) = Sk=1n[ C0kфk(j)]; j = 1…N (4.2.21)

Здесь S0(j) – идеальный (незашумленный) сигнал; фk(j), k = 1…n , - известные функции, составляющие идеальный сигнал; h(j) – случайный сигнал с известным законом распределения; а – коэффициент при случайном сигнале, с который условно назовем амплитудой помехи. В качестве базиса для разложения f(j) примем функции фk(j), входящие в идеальный сигнал (21). Целью разложения в этом случае является как можно более точное определение параметров C0k незашумленного сигнала, то – есть фильтрация помех. В результате разложения получим

S(j) = Sk=1n[Ck фk(j)] (4.222)

где коэффициенты разложения Ck вычислены по соотношениям (17) и (20).

Рис.3


Рис.4

Рис.4




















Рис.3
























Рис.3




На Рис. 4.3 представлено функциональное пространство H (для N = 3), плоскость Sn (для n = 2), а также все величины, входящие в (21) и (22). Вектор aSh(j) является проекцией вектора помехи ah(j) на плоскость Sn, то – есть разложением помехи по базису фk(j). Из треугольника а1а2а3 (рис.3) получим соотношение:

aSh(j) + d(j) = ah(j). (4.2.23)

Поскольку вектор ошибки d (отрезок a1a2)ортогонален плоскости Sn , то треугольник a1a2a3 является прямоугольным, поэтому, по аналогу теоремы Пифагора, из уравнения (23) можно получить:

a2ESh + Ed = a2Eh (4.2.24)

Здесь a2Ehэнергия помехи ah(j), a2ESh - энергия разложения aSh помехи ah(j), Edэнергия ошибки разложения d(j) = f(j) – S(j). Фильтрующие свойства разложения можно характеризовать коэффициентом фильтрации Кф, отражающим уменьшение модуля помехи при проектировании ее на базисную гиперплоскость:

Кф = [aSh(j)]/[ah(j)] = [Sh(j)]/[h(j)]. (4.2.25)

Из прямоугольного треугольника а1а2а3 (Рис.3) и (25) можно получить соотношение:

Кф = Rad {ESh/Eh } , (4.2.26)

где Eh – энергия случайного сигнала h(j), ESh – энергия разложения Sh(j) этого сигнала. При достаточно больших N эти величины обладают статистической устойчивостью и, согласно [Л.4.4], их осредненные значения могут быть оценены теоретически. По определению, энергия случайного сигнала h(j)

Eh = Sj=1 N h2(j)

Осредним эту величину по ансамблю реализаций h(j):

Eh = M{Sj=1 N h2(j)} = Sj=1 N M{h2 (j)} = Sj=1 N Dh = NDh = N(Gh)2 (4.2.27)

Здесь Dh и Gh - дисперсия и среднеквадратическое отклонение случайной величины h(j) соответственно. При известном (или предполагаемом) законе распределения эти величины известны. Так, для равномерного закона распределения с единичным диапазоном

Dh = 1/12 = 0.083; Gh = 1/2Rad{3} = 0.289 (4.2.28)

Разложение случайного сигнала h(j) представим:

Sh (j) = Sj=1 N Chk фк(j) (4.2.29)

Энергия этого разложения, согласно (12), имеет вид:

Esh = Sk=1 n Chk Vhk (4.2.30)

Величины Chk и Vhk связаны соотношениями (9) и (17). Введем матрицу Q, обратную U из (17), тогда

Q = U-1 = [Qji]; j = 1…n , i = 1…n . (4.2.31)

Решение уравнений (9) теперь можно записать:

Chk = Sr=1 n Qkr Vhr

Подставив это в (30), получим:

Esh = Sk=1 n Vhk * Sr=1 n Qkr Vhr = Sk=1 n Sr=1 n Qkr Vhr Vhk (4.2.32)

Осредним (32) по ансамблю реализаций:

Esh = M[Esh] = Sk=1 n Sr=1 n Qkr M[Vhr Vhk] (4.2.33)

В соответствии с (20) вычислим

M[Vhr Vhk] = M[Sj=1 N h(j) фk(j)* Si=1 N h(i) фr(i)] =

= Sj=1 N фk(j) Si=1 N фr(i) M[h(j)h(i)] (4.2.34)

Если случайные величины h(j) и h(i) независимы, то

M[h(j)h(i)] = { Dh …i = j

{ 0…i не =j (4.2.35)

Подставив (35) в (34) с учетом (20), получим:

M[Vhr Vhk] = Dh Sj+1 N фk(j) фr(j) = Dh Urk (4.2.36)

Подставим (36) в (33):

Esh = Dh Sk=1 n Sr=1 n Qkr Urk (4.2.37)

Внутренняя сумма в (37) является элементом произведения матриц Q и U, которые, согласно (31), взаимно обратны, поэтому

Sr=1 n Qkr Urk = Qkk Ukk = 1 (4.2.38)

Используя (38) , из (37) получим:

Esh = Dh n (4.2.39)

Подставив в выражение (26) для коэффициента фильтрации осредненные значения (27) и (39), получим:

Кф = Rad(n/N) (4.240)



4.2.5. Оценка уровня помех в обрабатываемой реализации

При обработке экспериментальных реализаций представляет интерес уровень помех, присутствующих в этих данных. Для сигналов, которые можно представить в виде зашумленных сигналов (21), можно оценить амплитуду помехи (а). Из соотношений (24) – (26) можно выразить амплитуду помехи, присутствующую в эксперименте:

a2 = Ed/(Eh – Esh ) = Ed/Eh(1 – Кф2)

Энергия ошибки разложения Ed вычисляется в процессе разложения обрабатываемой реализации. Величины Eh, Esh, и Кф оцениваются выражениями (27), (39) и (40). Таким образом:

a2 = Ed/Dh(N – n) (4.2.41)

Напомним, что здесь Dh – дисперсия шума, N – число отсчетов в обрабатываемой реализации, n – количество базисных функций.



4.2.6. Оценка погрешностей разложения

Для зашумленных сигналов типа (21) в качестве погрешности dS(j) разложения естественно принять разность между разложением S(j) зашумленного сигнала и идеальным (незашумленным) сигналом S0(j).

Учитывая выражения (21) и (22) можно записать:

dS(j) = S(j) – S0(j) = Sk=1 n {(Ck – C0kk(j)} = Sk=1 n {dCk фk(j)} (4.2.42)

Здесь dCk = Ck – C0k – погрешности разложения в области параметров сигнала.

Из рассмотрения треугольника a0a1a3 (рис.3) можно получить:

dS(j) = S(j) – S0(j) = аSh(j) (4.2.43)

Вектор аSh(j) был определен выше, как проекция помехи на базисную гиперплоскость Sn, амплитуда помехи а оценивается соотношением (41).

На Рис.4.4 представлен треугольник а0а1а3 из Рис.4.3, расположенный в базисной гиперплоскости Sn. При обработке экспериментальных данных вектор S(j) становится известным, вектора же S0(j) и аSh(j) – неизвестны. Однако длина L вектора aSh(j), характеризующего ошибку разложения, может быть оценена с помощью соотношений (39), (40),(41):

L = aRad{ESh} = Rad{Edn/(N-n)} = КфRad{Ed/(1 – Кф2 )} (4.2.44)



















Рис.4


На Рис.4.4 приведена окружность (в общем случае гиперокружность) с центром в точке а3 и радиусом L (31). Истинное положение точки а3 (конца вектора S0(j)) лежит на упомянутой окружности, что и является

оценкой погрешности разложения.

Сравнивая соотношения (42) и (43) можно получить:

Sk=1 n dCk фк(j) = a Sk=1 n Chk фк(j) (4.2.45)

Отсюда

dCk = aChk , k = 1…n. (4.2.46)

Уравнения (46) показывает, что статистические характеристики ошибок коэффициентов разложения Ck могут быть оценены по статистическим характеристикам коэффициентов Chk разложения шума.

Дисперсию Chk с учетом (31) можно записать:

D[Chk] = M[Chk 2] = M[Sr=1 n QrkVhr* Si=1 n QikVhi =

(4.2.47)

=Sr=1 n Qrk Si=1 n QikM[VhrVhi]

Подставив сюда (36), получим:


D[Chk] = Dh Sr=1 n Qrk Si=1 n Qik Uri (4.2.48)

Внутренняя сумма в (48) является элементом единичной матрицы в силу (31) , поэтому:

Si=1 n Qik Uri = { 1…r=k }

{ 0…r не = k}

Отсюда и из (48)

D[Chk] = Dh Qkk (4.2.49)

Таким образом, дисперсию ошибок вычисления коэффициентов разложения сигналов типа (21) можно выразить из (46) и (49) следующим образом:

D[dCk] = a2 Dh Qkk (4.2.50)

4.2.7. Численная оптимизация по параметрам базиса

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

f(j) = S0 (j) + ah(j) , S0 (j) = Sk=1 nC0kфk (j, q0i ) , (4.2.51)

i = 1…m , j = 1…N

В качестве базисной системы функций для разложения сигнала (51) естественно принять

ф1(j,qi), ф2(j,qi) …фm(j,qi) , i = 1…m. (4.2.52)

Здесь qi – неизвестные параметры базисных функций, m – число этих параметров. Задав произвольную совокупность параметров qi, i = 1…m, произведем разложение сигнала (51) по базису (52), используя соотношения (17) и (20). В результате получим набор коэффициентов разложения и энергию ошибки, которые являются функциями принятой совокупности параметров qi:

Ck = Ck(q1…qm) , Ed = Ed(q1…qm) (4.2.53)

Произведя численную оптимизацию Ed(q1…qm) по параметрам qi, i = 1…m, получим:

Ck опт = Ck опт(q1опт…qmопт), Edопт = Edопт(q1опт…qmопт) (4.2.54)

Оптимальные параметры из (54)

qiопт , i = 1…m, Ck опт, k = 1…n, (4.2.55)

являются оценками параметров q0i , i = 1…m , и С0k , k = 1…n незашумленного сигнала (51).



4.2.8 Оценка погрешностей разложения методом статистического моделирования

Погрешности оценок (55) параметров q0i и C0k незашумленного сигнала S0(j) (50) можно произвести с помощью статистического моделирования. Задав значения параметров q0i и C0k и амплитуду шума а, сформируем ансамбль зашумленных реализаций, отличающихся одна от другой реализациями шума h(j). Разложение каждой реализации из ансамбля даст свою оценку (55) параметров q0i и C0k , что позволяет вычислить погрешности всех параметров для данной реализации:

dqi = q0i – qi опт, i = 1…m, dCk = C0k – Ckопт, k = 1…n (4.2.56)

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

4.3. Идентификация динамического звена по экспериментальному частотному годографу

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

Экспериментальный годограф W e(i w) представим в виде известной совокупности N точек на комплексной плоскости с параметром wk:

W e(i w k) = W ek = N ek exp( iф ek) = D ek + i E ek , k = 1…N (4.3.1)

Здесь W ek – комплексная, N ek – амплитудная, ф ek – фазовая, D ek – вещественная, Eek – мнимая экспериментальные частотные характеристики. Эти характеристики связаны соотношениями:

(N ek)2 = (D ek)2 + (E ek)2; tg(ф ek) = E ek /D ek (4.3.2)

Будем считать, что характеристики приведены к нулевой частоте, т.е. на частоте w = 0

W e = 1; N e = 1; D e = 1; ф e = 0; Ee = 0 (4.3.3)

В качестве аппроксимирующей модели примем линейное динамическое звено с передаточной функцией

W(p) = (b m p m + b (m-1) p (m-1) …+b 1 p + 1)/(a n p n …a 1 p + 1) (4.3.4)

Здесь известными будем считать порядок числителя (m) и знаменателя (n), а неизвестными – значения коэффициентов a i (i = 1…n) и b j (j = 1…m). Количество неизвестных коэффициентов равно ( n + m ). Подставив в W(p) p = w k , получим выражения частотных характеристик модели для частот экспериментального годографа:

W k = [A k +i B k ]/[U k + i V k ] = D k + i E k = N k exp(iф k ) (4.3.5)

Здесь

A k = 1 – b 2 w k 2 + b 4 w k 4 …+ b m w k m

B k = b 1 w k – b 3 w k 3 …+ b (m-1) w k (m-1)

U k = 1 – a 2 w k 2 + a 4 w k 4 …+ a n w k n (4.3.6)

V k = a 1 w k – a 3 w k 3 …+ a (n-1) w k (n-1)

D k = (A k U k + B k V k)/(U k 2 + V k 2)

E k = (B k U k – A k V k)/ (U k 2 + V k 2) (4.3.7)

N k 2 = D k 2 + E k 2

tg(ф k ) = E k /D k

На рис.4.5 представлена комплексная плоскость, на которой кривой 1 изображен годограф объекта We(iw), на котором отмечены экспериментальные точки Wek , соответствующие частотам wk,,k =1…N. Кривой 2 изображен годограф W(iw) модели (4) при произвольных значениях ai и bj и отмечены точки Wk , соответствующие тем же частотам. Вектор dWk характеризует отличие годографов объекта и модели:

dW k = W ek – W k (4.3.8)

Выбор параметров модели будем производить из условий близости годографов объекта и модели.



4.3.2 Совпадение годографов в экспериментальных точках

Условием совпадения годографов в экспериментальных точках являются очевидные соотношения

dW k = W ek – W k = 0. (4.3.9)

Комплексное уравнение (9) эквивалентно двум вещественным, которые, с учетом (7), можно записать:

D ek = D k = (A k U k + B k V k)/(U k 2 + V k 2) (4.3.10)

E ek = E k = (B k U k – A k V k)/ (U k 2 + V k 2)

Поскольку одной экспериментальной точке соответствуют два уравнения, а количество искомых параметров передаточной функции модели (4) равно (n + m), то для принципиальной разрешимости уравнений (10) количество N экспериментальных точек должно быть:

k = 1…N ; N = (n + m)/2 (4.3.11)

Из (10) и (6) видно, что система уравнений (10) нелинейна, что затрудняет ее решение. Однако легко показать, что можно построить эквивалентную линейную систему уравнений. Для этого представим уравнение (9) с учетом (1) и (5) в следующем виде:



dW k = D ek + i E ek – (A k +i B k )/(U k + i V k ) = 0 (4.3.12)

Умножив (12) на отличный от нуля множитель (U k + i V k) и разделив вещественную и мнимую части, получим линейную систему уравнений, эквивалентную системе (10):

D ek U k – E ek V k – A k = 0 ; E ek U k + D ek V k – B k = 0 ; (4.3.13)

k = 1…(n+m)/2

4.3.3. Оптимизация модели при избытке экспериментальных точек

Если количество N экспериментальных точек в годографе превышает величину (11), то провести годограф модели через все эти точки невозможно. В этом случае необходимо ввести приемлемый критерий, характеризующий близость двух годографов на всей совокупности экспериментальных точек, и минимизировать его по параметрам модели. В качестве критерия близости можно было бы выбрать сумму квадратов модулей расхождений dW k (8):

I1=Sk=1 N[dW k] 2=S k=1 N(N ek 2 + D k 2 + E k 2 – 2D ek D k – 2E ek E k) (4.3.14)

Минимизация I1 приводит, как и в случае (10), к нелинейной системе уравнений для определения параметров модели. Попытаемся несколько видоизменить критерий близости годографов с целью получения линейной системы уравнений при его оптимизации. Из первого уравнения (12)

`dW k = D ek + i E ek – (A k +i B k )/(U k + i V k ) (4.3.15)

Умножим уравнение (15) на отличный от нуля комплексный множитель (U k + i V k ):

dH k = dW k (U k + i V k ) = (D ek U k – E ek V k – A k) + …

+ i (E ek U k + D ek V k – B k ) (4.3.16)

В качестве отличия годографов в экспериментальной точке примем dHk , а в качестве критерия близости I – сумму квадратов модулей:

I = S k=1 N[ dH k] 2 = S k=1 N{(A k 2 + B k 2 ) + N ek 2 (U k 2 + V k 2 ) – …

- 2D ek (A k U k + B k V k ) – 2E ek (A k V k – B k U k )} (4.3.17)

Отметим, что критерии близости I1 (14) и I (17) в нуль обращаются одновременно. Квадратичная форма I (17) является функцией коэффициентов a i , b j модели. Для минимизации ее приравняем нулю частные производные от I по этим коэффициентам:

dI/da i = 0; i = 1…n . dI/db j = 0 ; j = 1…m. (4.3.18)

В выражениях (6) параметры зависят только от четных или нечетных коэффициентов модели:

A k = 1 + S i=1 m/2(-1) i w k 2i b 2i , B k = S i=1 m/2(-1) i+1 w k (2i-1) b (2i-1) ,

(4.3.19)

U k = 1 + S j=1 n/2(-1) j w k 2j b 2j , V k = S j=1 n/2(-1) j+1 w k (2j-1) b (2j-1) .

Выпишем отличные от нуля производные, которые потребуются для (18):

dA k /db 2q = w 2q (-1) q , dB/db (2q-1) = w (2q-1) (-1) (q-1) , q = 1…m/2(4.3.20)

dU k /da 2r = w 2r (-1) r , dV k/da (2r-1) = w (2r-1) (-1) (r+1) , r = 1…n/2

Уравнения (18) с учетом (20) можно записать:

dI/da 2q = 2(-1) q S k=1 N(N ek 2 U k - D ek A k + E ek B k )w k 2q = 0

dI/da 2q-1 = 2(-1) q+1 S k=1 N(N ek 2 V k - D ek B k E ek A k )w k (2q-1) =0 (4.3.21)

dI/db 2r = 2(-1) r S k=1 N(A k - D ek U k - E ek V k )w k 2r = 0

dI/db 2r-1 = 2(-1) r+1 S k=1 N(B k - D ek V k + E ek U k )w k (2r-1) = 0

q = 1…m/2 , r = 1…n/2 .

Система уравнений (21) является линейной относительно параметров модели (4). Парамнтры ai и bj , вычисленные из (21), минимизируют выбранный критерий I (17) близости годографов, и , таким образом, идентифицируют испытуемое динамическое звено.





4.4. Идентификация динамического звена по экспериментальному переходному процессу

4.4.1. Идентификация аналогового динамического звена

Подавая на вход испытуемого динамического звена ступенчатое возмущение, получим экспериментальный переходный процесс Хэh) ввиде совокупности N его отсчетов, взятых с шагом h:


Хэ(t) = Xэ(k) , t = kh, k = 1…N (4.4.1)


В качестве аппроксимирующей модели примем линейное динамическое звено с передаточной функцией


X(p)/Y(p) = W(p) = [bm pm +… b1 p + b0 ] / [pn + an-1 pn-1 +…a1 p + a0]


При ступенчатом входном воздействии Y(p) = 1/p и нулевых начальных условиях переходный процесс описывается соотношением [Л.4.5] :


X(t) = D0 + Sj=1 n [ Dj exp(pj t) ] (4.4.2)


Здесь pj – корни характеристического уравнения


A(p) = pn + an-1 pn-1 +… a1 p +a0 = 0



Выражение (2) показывает, что в качестве базиса для разложения переходного процесса (1) следует принять:


ф0(k) = 1, ф1(k) = exp(s1hk), … фn(k) = exp(snhk) (4.4.3)


Представим


Xэ(k) = S(k) + d(k) = Sj=1n [Cj exp(sjhk)] +d(k)


Здесь S(k) – разложение, d(k) – ошибка разложения.

Минимизация энергии ошибки Ed = Sk=1N[d(k)2] по коэффициентам Cj при заданных значениях парметров sj базиса (3) определяется, как показано в разделе (4.2.1), системой уравнений (4.2.9)

с обозначениями (4.2.20):


Sjn[Cj Uij ] = Vi , i = 1…n

Uij = Sk=1Ni(k)*фj(k)] , Vi = Sk=1Ni(k)*Xэ(k)]


Коэффициенты разложения Cj и энергия ошибки Ed являются функциями от параметров si базиса (3) :


Cj = Cj(s1 …sn) ; Ed = Ed (s1… sn )


Численная оптимизация Ed по параметрам sj позволят найти набор параметров siопт и коэффициентов Cj опт который позволяет сформировать оптимальную сумму


Sопт(t) = C0опт + Sj=1n[Cjопт exp(sjопт t)] , (4.4.4.)


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

выражение (2). Применив к (4) операцию Лапласа, получим изображение оптимальной суммы:


Sопт(p) = C0опт /p + Sj=1n[Cjопт /(psjопт )]


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


Sопт(p) = [bmопт pm +…b1опт p + b0опт ] / [ pn + …a1опт p + a0опт] p . (4.4.4)


Коэффициенты biопт и ajопт из (4) идентифицируют испытуемую систему





4.4.2. Идентификация дискретного динамического звена


Пусть испытуемый дискретный динамический объект описывается разностным уравнением [Л.4.6.]:


X(k+n) + Si=0n-1ai x(k+i) = Si=0mbi y(k+i) (4.4.5)


Здесь y(k) и x(k) – входная и выходная дискретные последовательности соответственно.

Порядок разностей n и m <= n левой и правой частей уравнения (5) будем считать известными, а коэффициенты ai и bj – неизвестными, подлежащими определению из результатов испытаний. Испытания заключаются в подаче на вход объекта некоторой дискретной последовательности y(k) и экспериментальном определении выходной дискретной последовательности x(k). Обе последовательности фиксируются (запоминаются) на некотором интервале 0 <= k <= N, как представлено на рис.4.6. На рисунке выделена произвольная группа отсчетов длиной (n+1), входящая в уравнение (5), которое должно удовлетворяться для любой группы отсчетов соответствующей длины. В уравнении (5) n+m+1 неизвестных коэффициентов ai , bj , поэтому для определения всех коэффициентов необходимо сформировать столько – же уравнений. Для этого необходимое наименьшее количество отсчетов в экспеиментальной реализации


Nmin = 2n + m +1 (4.4.6)


Если N > Nmin , то удовлетворить все уравнения типа (5) , которые можно сформировать на такой реализации, невозможно. Для подобного случая сформируем невязку:


L(k) = x(k+n) + Si=0n-1 ai x(k+i) - Si=0m bi x(k+i) (4.4.7)


В качестве критерия, характеризующего удаленность уравнения (5) от экспериментальной совокупности отсчетов, примем энергию

невязки:

EJ = Sk=0N-n L2(k) (4.4.8)


Наилучшей совокупностью кэффициентов ai , bj будем считать ту, которая минимизирует EJ , поэтому, из (8):


DEJ/Dai = Sk=0N-n 2L(k)DLk/Dai = 0 , i = 0…n-1

(4.4.9)

DEJ/Dbj = Sk=0N-n 2L(k)DLk/Dbj = 0 , j = 0…m


Из (7):

DLk/Dai = x(k+i), DLk/Dbj = x(k+j) (4.4.10)


Подставив (10) в (9) и введя обозначения:


qij = Sk=0N-n x(k+i)x(k+j) = qij


pij = Sk=0N-n y(k+i)y(k+j) = pij (4.4.11)


rij = Sk=0N-n x(k+i)y(k+j) = rij ,


получим:


Si=0n-1 aiqij – Si=0m birji = - qnj , j = 0…n-1

(4.4.12)


Si=0n-1 airij – Si=0m bipij = - qnj , j = 0…n-1


Значения ai , bj , вычисленные из (11) и (12) , минимизируют невязку (8), то-есть идентифицируют испытуемую дискрктную систему.


4.5. Упражнения


1) Сформировать ансамбль модельных реализаций зашумленного гармонического сигнала, как в разделе 4.2.4.:

fm(j) = S0(j) + ahm(j) на интервале 0<= j <= N

Здесь S0(j) = cos(bj) - идеальный сигнал

hm(j) – m – тая реализация случайной помехи с равномерным законом распределения

1<= m <= M , M – число реализаций в ансамбле

b = 2 pi r(1 + u)/N - частота идеального сигнала

r – целое число периодов в модельной реализации

0 <= u <= 1 , u – дробная часть периода в модельной реализации

Номинальные значения существенных параметров реализаций :

N = 100 ; r = 2 ; u = 0.3 ; a = 1 ; M = 20 .

Разложить реализации ансамбля с номинальными значениями параметров по базису:

ф0(j) = 1, ф1(j) = cos(dj), ф2(j) = sin(dj).

При численной оптимизации по частоте базиса d принять диапозон поиска dopt : 0.95b <= d < = 1.06b .

Допустимую погрешность по dopt принять равной 10-5 b .

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



3) Повторить пункт 2) , проварьировав амплитуду помехи в диапозоне 0.2 <= a <= 2 .

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



4) Повторить пункт 2), проварьировав число отсчетов в реализации в диапозоне: 10 <= N <= 200

Построить зависимости среднеквадратических значений погрешностей от числа отсчетов.



5) Повторить пункт 2) , проварьировав число целых периодов в реализации в диапозоне 10 <= r <= 200 .

Построить зависимости среднеквадратических значений погрешностей от r.



6) Повторить пункт 2) , проварьировав дробную часть периода в реализации в диапозоне 0 <= u <= 0.5.

Построить зависимости среднеквадратических значений погрешностей от u.



7) Повторить пункт 2) , проварьировав число реализаций в ансамбле в диапозоне 10<= M <= 100 .

Построить зависимости среднеквадратических значений погрешностей от числа реализаций M.



Литература

[Л.4.1] Эйкхофф П. Основы идентификации систем управления. – М.: Мир, 1951.- 683 стр.

[Л.4.2] Трахтман А.М. Введение в обобщенную теорию сигналов. – М.: Советское радио, 1972. - 352 стр.

[Л.4.3] Крылов А.Н. Лекции о приближенных вычислениях. – М-Л.: ГИТТЛ, 1950.- 398 стр.

[Л.4.4] Вентцель Е.С. Теория вероятностей. – М.: ФМ, 1962.-564 стр.

[Л.4.5] Диткин В.А. , Кузнецов П.И. Справочник по операционному исчислению. – М-Л.: ГИТТЛ, 1951.- 255стр.

[Л.4.6] Кузовков





29