авторефераты диссертаций БЕСПЛАТНАЯ БИБЛИОТЕКА РОССИИ

КОНФЕРЕНЦИИ, КНИГИ, ПОСОБИЯ, НАУЧНЫЕ ИЗДАНИЯ

<< ГЛАВНАЯ
АГРОИНЖЕНЕРИЯ
АСТРОНОМИЯ
БЕЗОПАСНОСТЬ
БИОЛОГИЯ
ЗЕМЛЯ
ИНФОРМАТИКА
ИСКУССТВОВЕДЕНИЕ
ИСТОРИЯ
КУЛЬТУРОЛОГИЯ
МАШИНОСТРОЕНИЕ
МЕДИЦИНА
МЕТАЛЛУРГИЯ
МЕХАНИКА
ПЕДАГОГИКА
ПОЛИТИКА
ПРИБОРОСТРОЕНИЕ
ПРОДОВОЛЬСТВИЕ
ПСИХОЛОГИЯ
РАДИОТЕХНИКА
СЕЛЬСКОЕ ХОЗЯЙСТВО
СОЦИОЛОГИЯ
СТРОИТЕЛЬСТВО
ТЕХНИЧЕСКИЕ НАУКИ
ТРАНСПОРТ
ФАРМАЦЕВТИКА
ФИЗИКА
ФИЗИОЛОГИЯ
ФИЛОЛОГИЯ
ФИЛОСОФИЯ
ХИМИЯ
ЭКОНОМИКА
ЭЛЕКТРОТЕХНИКА
ЭНЕРГЕТИКА
ЮРИСПРУДЕНЦИЯ
ЯЗЫКОЗНАНИЕ
РАЗНОЕ
КОНТАКТЫ


Pages:     | 1 || 3 |

«ФЕДЕРАЛЬНОЕ АГЕНТСТВО ПО ОБРАЗОВАНИЮ ГОСУДАРСТВЕННОЕ ОБРАЗОВАТЕЛЬНОЕ УЧРЕЖДЕНИЕ ВЫСШЕГО ПРОФЕССИОНАЛЬНОГО ОБРАЗОВАНИЯ «САМАРСКИЙ ...»

-- [ Страница 2 ] --

Приведенный алгоритм усреднения реализует сглаживание про цесса линейным фильтром с длиной памяти T = m t и амплитудно фазовой частотной характеристикой вида:

mt Wф ( j ) = (3.20) sin.

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

Следует заметить [64], что получаемые при таких процедурах ус ) реднения оценки y (t ) всегда являются смещенными, т.к. любая про цедура сглаживания соответствует прохождению зашумленного сиг нала через некоторый фильтр, отделяющий низкочастотный полез ный сигнал от более высокочастотной помехи.

Рисунок 3. Частотная характеристика фильтра Wф ( j ) и частотные спектры полезно го сигнала S п с ( ) и помехи S п ( ) На рисунке 3.7. приведена графическая интерпретация процесса фильтрации полезного сигнала со спектром S п с ( ) и высокочастот ной помехи со спектром S п ( ) низкочастотным фильтром с характе ристикой Wф ( j ). При малой полосе пропускания фильтра помеха лучше отфильтровывается, но подавляется и часть полезного сигнала.

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

Также часто используются и другие методы сглаживания – метод сглаживания четвертыми разностями, метод сглаживания с использо ванием разложения в ряд Фурье, разложения с использованием поли номов Чебышева [8, 64]. Применение рядов Фурье и полиномов Че бышева дает лучшее качество, но и требует значительно больших объёмов вычислений.

3.5 Идентификация объектов с помощью частотных характери стик Частотный метод идентификации линейных систем основан на работах Найквиста и Боде и использует в качестве исходной частот ные или спектральные характеристики.

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

Для линейного стационарного объекта вход-выходное соотноше ние определяется через частотную передаточную функцию:

y ( j ) = W ( j ) u ( j ), (3.21) где u ( j ) = F {u (t )} = j t u (t ) e dt - частотный спектр (преобразо вание Фурье) входного сигнала объекта;

y ( j ) = F {y (t )} - частотный спектр выходного сигнала;

W ( j ) = F {w (t )} - частотная передаточ ная функция объекта.

Из (3.21) следует, что частотную характеристику объекта y ( j ) W ( j ) = можно найти экспериментальным путем на основе u ( j ) частотных спектров измеренных входных и выходных сигналов u (t ) и y (t ). Основной сложностью при таком подходе является невозмож ность формирования входного воздействия u (t ), частотный спектр которого был бы непрерывным на бесконечном интервале изменения частот. Частотная характеристика объекта W ( j ) также должна быть непрерывной во всей полосе частот, что обеспечить технически практически невозможно. Поэтому, вследствие наблюдения сигналов u (t ) и y (t ) только на ограниченном отрезке времени, возникают зна чительные ошибки их измерения.

В соответствии с этим, в реальных условиях воздействие поли гармонических сигналов в широком диапазоне изменения частот за меняют последовательным применением моногармонических воздей ствий u ( t ) = u 0 sin i t с разными частотами = i и исследуют реак цию на них. На выходе объекта в установившемся состоянии будут наблюдаться гармонические колебания той же частоты y ( t ) = y m ( i ) sin[ i t + ( i )]. (3.22) В этом случае частотная передаточная функция, представленная в комплексном виде, определяется зависимостью:

y ( j ) y m ( i ) e j ( i ) = W ( j ) e j ( i ), (3.23) W ( j ) = = u ( j ) u где W ( j ) = mod W ( j );

( ) = arg W ( j ) - амплитудно- и фазо частотные характеристики объекта. В соответствии с (3.23) модуль функции W ( j i ) при заданной частоте тестового сигнала i вычис ляется по формуле:

y m ( i ) W ( j i ) =, (3.24) u а аргумент W ( j i ) определяется величиной фазового сдвига выход ных колебаний исследуемого объекта.

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

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

Частотная характеристика объекта представляется в следующем виде:

W ( j ) = W ( j ) e j ( ) = W ( j ) cos ( ) + + j W ( j ) sin t = c ( ) + jd ( ), (3.25) где c ( ) = W ( j ) cos ( ), d ( ) = W ( j ) sin ( ) являются под лежащими определению параметрами частотных характеристик объ екта W ( j ) = {c 2 ( ) + d 2 ( )} и ( ) = arctg d ( ) для ряда час c ( ) тот i.

Параметры c(i ) и d ( i ) можно найти, подавая на вход объекта тестовый сигнал u (t ) = u 0 sin i t и применяя фильтрацию Фурье к выходному сигналу y (t ). Такая процедура фильтрации реализуется умножением y (t ), соответственно, на sin i t и cos i t и усреднени ем по целому k числу периодов:

T c ( i ) = y (t ) sin i tdt, u 0T T d ( i ) = y (t ) cos i tdt, (3.26) u 0T где T = kTi, Ti =.

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

Когда выход объекта y (t ) искажается аддитивным белым шумом (t ), M [ ] = 0, то в результате фильтрации вместо детерминирован ных величин c(i ) и d (i ) будут определяться случайные величины ~ ~ c ( i ) и d ( i ) [74]:

T ~ (, T ) = c ( ) + 2 (t ) sin tdt ;

сi i i u 0T (3.27) T ~ d ( i, T ) = d ( i ) + (t ) cos i tdt.

u 0T Определение математического ожидания по всему ансамблю на блюдений показывает, что оценки являются несмещенными:

T ~ ] = M [c ( ) ] + 2 M [ (t ) ]sin tdt = c ( );

M [c i i i u 0T (3.28) [] T ~ M d = M [d ( i ) ] + M [ (t ) ]cos i tdt = d (i ).

u 0T Оценим дисперсию оценок на частоте i :

~ 2 (c ) = M [( c c ) 2 ] = M [( (t ) sin i tdt ) 2 ] = 2 uT (3.29) 4 TT = 2 2 M [ (t ) ( ) ]sin i t sin idtd u0 T 0 В (3.29) M [ (t ) ( ) ] - корреляционная функция аддитивного шума, и ее величина которой равна:

2, t = 0;

K (t ) = M [ (t ) ( ) ] = (3.30) 0.

С учетом соотношения (3.30), вычисляя двойной интеграл в (3.29), можно определить дисперсию оценок c(i ) :

4 T 2 T (c ) = 2 2 sin i t sin idtd = u 2T 2 2 = u 2T, 2 (3.31) u0 T 0 0 где - дисперсия шума;

u 0 - амплитуда входного сигнала [74].

Аналогично определяется дисперсия оценок d (i ) :

(d ) = (3.32).

u0 T Полученные результаты (3.31) и (3.32) характеризуют зависи мость дисперсии оценок параметров частотных характеристик иден тифицируемого объекта от величины дисперсии аддитивного шума, интенсивности входного сигнала и величины выборки и показывают малую чувствительность алгоритмов гармонического анализа сигна лов к содержащимся в контуре идентификации помехам типа «белого шума».

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

3.6 Корреляционные методы Корреляционные методы являются основными в случае исполь зования непрерывных (невыборочных) сигналов. Они широко ис пользуются для решения задач идентификации систем управления, обнаружения сигналов, технических измерений, реализации иссле дуемых случайных сигналов и т.д. Корреляционные методы основы ваются на дифференциальной аппроксимации и применяют опреде ленный оператор (в частности, умножение на функцию времени и ин тегрирование произведения) к входным и выходным сигналам [19, 20, 53, 54, 64, 65, 74].

Рассмотрим линейный стационарный одномерный объект, на ко торый воздействуют непрерывные стационарные сигналы u (t ) и (t ), и задача идентификации состоит в оценке весовой функции объекта w(t ) по наблюдениям y (t ) и u (t ) на некотором интервале времени [0, T ].

Рисунок 3. Схема идентификационного эксперимента по оценке импульсной весовой функции объекта Так как объект линейный, то для отдельных реализаций случай ных процессов выходной сигнал описывается интегралом свертки:

y ( t ) = w ( ) u ( t ) d + ( t ) = w ( t ) * u ( t ) + ( t ). (3.33) Выражение (3.33) связывает единичные реализации случайных вход ного и выходного процессов.

Для перехода к детерминированным величинам проведем умно жение обеих частей (3.33) на u (t ) и применим операцию матема тического ожидания, в результате чего получим выражение:

M [u ( t ) y (t ) ] = w ( ) M [u ( t ) u (t ) ]d + + M [u (t ) ( t ) ]. (3.34) Учитывая, что соотношения K uu ( ) = M [u (t )u (t ) ];

K uy ( ) = M [u (t ) y (t ) ] (3.35) определяют соответственно автокорреляционную функцию входного сигнала и взаимную корреляционную функцию, перейдем к уравне нию относительно корреляционных функций сигналов:

K uy ( ) = w ( ) K uu ( ) d + 0 = w ( ) K uu ( ) + 0. (3.36) Если выполняются условия физической реализуемости системы, т.е.

w(t ) = 0 при t 0, и сигналы u(t) и (t ) не коррелируемы, то (3.36) трансформируется в уравнение, называемое уравнением Винера – Хопфа:

K uy ( ) = w ( ) K uu ( ) d. (3.37) Уравнение (3.37) связывает детерминированные величины – корреля ционные функции сигналов, - и его решение позволяет оценить им пульсную весовую функцию линейной стационарной системы по критерию минимума среднеквадратической ошибки:

1T 1T J = [ y (t ) y M (t )] dt = [ y (t ) w ( )u (t ) d ]2 dt. (3.38) T0 T0 Решение уравнения Винера-Хопфа (3.37) аналитическими мето дами вызывает значительные сложности, связанные, во-первых, с трудностью получения аналитического описания корреляционных функций по реализациям входо-выходных сигналов, полученных экс периментально в режиме нормальной эксплуатации объекта. Во вторых, теоретически корреляционные функции должны определять ся по бесконечным выборкам. На практике же они оцениваются на конечном интервале времени [0, T ], что приводит к приближенным оценкам истинных корреляционных функций. Появление случайных ошибок в оценках корреляционных функций влечет существенные погрешности при нахождении w(t ), и задача решения уравнения Ви нера-Хопфа аналитическими методами становится некорректной [65].

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

Наиболее простой способ идентификации w(t ), позволяющий обойтись без решения уравнения свертки (3.37), заключается в ис пользовании в качестве входного воздействия сигнала, корреляцион ная функция которого близка по свойствам к импульсной функции Дирака K uu ( ) = c ( ), т.е. сигнала типа белого шума, где c – интен сивность шума. В этом случае ядро интеграла Винера-Хопфа и дает приближенную оценку весовой функции объекта K uy (t ) cw (t ). (3.39) Отсюда по значению взаимной корреляционной функции и определя ется w(t ) объекта:

w (t ) K uy (t ). (3.40) c Таким образом, подавая на вход линейного стационарного объек та сигнал типа белого шума и находя его взаимную с выходным сиг налом корреляционную функцию, можно непосредственно получить идентифицируемую импульсную переходную характеристику (ИПФ).

Пример 3. Проиллюстрируем определение импульсной весовой функции апериодического объекта первого порядка с запаздыванием, переда 10 10 p точная функция которого равна W ( p) =, по результатам e 3p + проведения эксперимента при входном сигнале типа белого шума.

s1=tf([10],[3 1], 'td', 10) % непрерывная передаточная функция объекта T_izm=60;

% интервал измерений T_end=30;

% интервал оценивания dt=0.05;

% шаг дискретизации t_izm=0:dt:T_izm;

% массив дискретного времени измерений t=0:dt:T_end;

% массив дискретного времени оценивания N_izm=length(t_izm);

% размер выборки интервала измерений N=length(t);

% размер выборки интервала оценивания u=randn(N_izm,1);

% массив значений входного воздействия y=lsim(s1,u,t_izm);

%массив значений выходного воздействия Ruu=xcorr(u,u,'biased');

% вычисление корреляционной функции на интер вале измерений по одной реализации Ruy=xcorr(y,u,'biased');

% вычисление взаимной корреляционной функции по одной реализации % вычисление усредненных авто- и взаимной корреляционных функций по 10 реализациям Ruu_mid=0;

% начальное значение автокорреляционной ф-ции Ruy_mid=0;

% начальное значение взаимной корреляционной ф-ции for k=1: u=randn(N_izm,1);

% массив значений входного воздействия y=lsim(s1,u,t_izm);

%массив значений выходного воздействия Ruu_mid= Ruu_mid+xcorr(u,u,'biased');

Ruy_mid= Ruy_mid+xcorr(y,u,'biased');

end;

Ruu_mid= Ruu_mid/10;

% усредненное значение автокоррел. ф-ции Ruy_mid= Ruy_mid/10;

% усредненное значение взаимной коррел. ф-ции tau=- T_end:dt: T_end;

figure (1);

plot (tau, Ruu(N_izm-N+1: N_izm +N-1));

figure (2);

plot (tau, Ruu_mid(N_izm-N+1: N_izm +N-1));

w= impulse(s1,t);

figure (3);

plot (t, w, t, 1/dt* Ruy (N_izm: N_izm +N-1));

figure (4);

plot (t, w, t, 1/dt* Ruy_mid (N_izm: N_izm +N-1));

На рисунках 3.9, 3.10 графически представлены результаты расчета.

1. 1. 0. 0. Ruu mid(t) 0. 0. Ruu(t) 0. 0. 0. 0. -0. -0. -30 -20 -10 0 10 20 -30 -20 -10 0 10 20 Время, с а б Рисунок 3. Корреляционная функция входного сигнала типа белый шум, получен ная по одной реализации (а) и усредненная по десяти реализациям (б) 3. 3. Импульсная весовая характеристика Импульсная весовая характеристика 3 2. 2. w(t), Ruy mid(t) w(t), Ruy(t) 1. 1. Взаимная Усредненная взаимная корреляционная функция корреляционная функция 0. 0. -0. -0. - 0 5 10 15 20 25 0 5 10 15 20 25 Время, с Время, с а б Рисунок 3. Истинная (аналитическая) импульсная переходная функция объекта и вза имная корреляционная функция входного и выходного сигналов, полу ченная по одной реализации (а) и усредненная по десяти реализациям (б) Существование определенных погрешностей между полученны ми взаимными корреляционными функциями и аналитической весо вой характеристикой объясняется, во-первых, невозможностью сформировать на входе объекта идеальный белый шум, а во-вторых, использованием не истинных корреляционных функций процесса, а их оценок, полученных на ограниченных временных интервалах.

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

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

Если в качестве входного сигнала используется сигнал, отличный от белого шума, то для определения оценки ИПФ используется под ход, заключающийся в численном решении дискретного аналога уравнения Винера-Хопфа (3.37) по оценкам корреляционных функ ций K uy ( ) и K uu ( ). В этом случае, в исходном уравнении (3.37) бесконечный предел интегрирования заменяется на конечный T, ко торый разбивается на N частей с интервалом дискретизации t, ин теграл заменяется на конечную сумму:

N w(it ) K uu [( j i)t ]t, j = 0,1,... N 1.

K uy ( jt ) = (3.41) i = Отсюда получается система N алгебраических уравнений, которая в развернутом виде принимает вид:

... K uu (( N 1)t ) K uy (t ) K uy (0) K uu (0) K (t ) K (t ) K uu (( N 2)t ) K uu (0) = t uy uu......

K uu (( N 1)t ) K uu (( N 2)t ) K uy (( N 1)t ) K uu (0) w(0) w(t ) w(( N 1)t ) или K uy = t K uu w. (3.42) Для решения системы уравнений (3.42) приходится использовать специальные методы параметрической идентификации (будут рас смотрены далее), например, Гаусса, коллокации, или наиболее рас пространенный – метод наименьших квадратов, в соответствии с ко 1 торым w = K uu K uy.

t В (3.42) матрица K uu представляет собой квадратную симмет ричную матрицу из значений автокорреляционной функции входного сигнала. Она обладает плохой обусловленностью, вследствие чего малые изменения коэффициентов матрицы ведут к большим погреш ностям в решении, и полученное решение w(t ) часто далеко от ис тинного. В этом случае для получения более точного решения реко мендуется пользоваться методом аппроксимирующих функций с ис пользованием метода коллокации [29] или специальными методами регуляризации [44, 64, 66].

4 Методы параметрической идентификации Перемены в прогрессивной стране неизбеж ны. Перемены - это постоянная величина Бенджамин Дизраэли 4.1 Общий подход к оцениванию параметров Методы параметрической идентификации широко применяются при решении практических задач идентификации линейных и нели нейных систем. Вопросу оценивания параметров посвящено много работ, начиная с разработки К. Гауссом метода наименьших квадра тов в 1795г., и на данный момент существует множество различных подходов, способов и методов параметрической идентификации [13, 19, 62, 71, 74]. Вместе с тем, базовым подходом к определению пара метров моделей остается метод наименьших квадратов (МНК) [23, 36, 38, 67, 71]. Дальнейшее развитие использования МНК привело к по явлению различных модификаций этого метода: обобщенный МНК (марковские оценки), метод взвешенных наименьших квадратов, ме тод штрафных функций и др. Наряду с МНК, существуют другие ме тоды идентификации, находящие свое применение, например, методы на основе вероятностного подхода: оценки максимального правдопо добия, максимума апостериорной информации, байесовские оценки, алгоритмы стохастической аппроксимации и прочие [37, 48, 57, 58, 59, 61, 68, 74].

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

Общие требования к методам идентификации заключаются в том, что определяемые оценки параметров модели должны быть точными и достигаться достаточно быстро [74]. В соответствии с этим, методы оценивания параметров должны быть:

- формализуемыми в достаточно общем виде;

- легко реализуемыми и обеспечивающими приемлемую скорость сходимости;

- обеспечивающими получение оптимальных оценок искомых па раметров.

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

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

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

математическим ожиданием M [ ], смещением {M [ ] b}, ковариа цией cov[ ] = M [{ M [ ]}{ M [ ]}T ]. Ковариация представляет собой матрицу, главная диагональ которой состоит из оценок диспер сий оцениваемых параметров, а недиагональные элементы соответст вуют оценкам взаимных дисперсий параметров.

В статистическом плане оценки вектора параметров b должны обладать определёнными свойствами [1, 74]: быть несмещенными, состоятельными, эффективными и достаточными, т.е.:

1. оценка будет несмещенной, если для любого N математиче ское ожидание оценки каждого параметра равно его истинному зна чению M [ ] = b ;

2. оценка является состоятельной, если для сколь угодно малого 0 с ростом N практически наверняка значение близко к b. На дежность оценки при увеличении выборки растет:

0 lim P[ b ] = 0 ;

N 3. оценка будет эффективной, если имеет наименьшую дис персию по сравнению с любыми другими оценками данного пара метра:

cov[ ] = M [{ b}{ b}T ] M [{ b}{ b}T ] = cov[ ].

4. оценка является достаточной, если для всех остальных оценок условная плотность вероятности P( / ) не зависит от b.

Для того, чтобы полученные оценки обладали желательными свойствами, делаются следующие допущения о случайных ошибках :

1. ошибки i, i = 1,2...N являются случайными величинами;

2. математическое ожидание ошибок i, i = 1,2... N равно нулю:

M [ ] = 0 ;

3. дисперсия ошибок i, i = 1,2...N постоянна D[ i ] = 2 ;

4. различные значения i независимы друг от друга:

0, при i j cov( i, j ) =, при i = j.

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

Пусть выходной сигнал объекта у состоит из аддитивной смеси отклика на входное воздействие u и шума :

y ( j ) = y (u,, b, j ) j = 1,2,...N. (4.1) Здесь b = [b0, b1,...bm ]T – вектор параметров объекта размерностью [(m + 1) 1];

j – номер наблюдаемого измерения выходной величины у;

N – размер выборки измерений ( N m + 1).

Будем считать, что сигналы u и у могут быть измерены точно, помеха характеризуется нулевым математическим ожиданием M [ ] = 0 и ковариационной матрицей H:

M [ (1) (1)]... M [ (k ) (1)] [ ] = = H.

cov[ ] = M T...

M [ (1) (k )]... M [ (k ) (k )] (4.2) Пусть в результате эксперимента имеется выборка измерений u = [u (1), u (2).......u ( N )] T входного и выходного сигналов y = [ y (1), y (2)....... y ( N )]. Для получения на основе эксперименталь T ных данных числовых оценок = {u (1), u (2),..., u ( N ), y (1), y (2),..., y ( N )} искомых параметров ис пользуется модель, аналогичная (4.1):

y M ( j ) = y M (u,0,, j ) j = 1,2,...N, (4.3) где = [ 0, 1,... m ] - вектор параметров модели размерностью T [(m + 1) 1].

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

m y M ( j) = i ui ( j) j = 1,2,...N, (4.4) i = которая в векторной форме принимает вид y M = U, (4.5) где матрица U представляет совокупность значений входных воздей ствий на объект:

u0 (1)... um (1) U =....

u0 ( N ) um ( N ) 4.2 Оценивание параметров объектов по методу наименьших квадратов.

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

по имеющимся выборочным данным наблюдений за входным и вы ходным сигналами с интервалом дискретизации t требуется оце нить значения параметров, обеспечивающих минимум величины функционала невязки между модельными и фактическими данными N J = ( y U ) ( y U ) = e e = e 2 ( j ).

T T (4.6) j = Здесь величина e( j ) = y ( j ) y M ( j ), j = 1,2...N, (4.7) представляет невязку, определённую как разность между выходом исследуемого объекта y = [ y (t ), y (2t ),.. y ( Nt )]T и реакцией, вы численной по математической или физической модели объекта y M = [ yM (t ), y M (2t ),..., yM ( Nt )]T. Невязка складывается из неточ ностей структуры модели, погрешностей измерений и неучтённых взаимодействий среды и объекта. Однако, независимо от происхож дения возникающих ошибок, МНК минимизирует сумму квадратич ной невязки для дискретных значений.

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

Оценка по МНК *, минимизирующая критерий (4.6), находится из условия существования минимума функционала:

J = min J = J = *. (4.8) Важным свойством оценок по МНК является существование только одного локального минимума, совпадающего с глобальным. Поэтому оценка * является единственной. Ее значение определяется из усло вия экстремума функционала (4.6):

J = 2U T ( y U * ) = 0, (4.9) = * откуда следует соотношение, определяемое систему нормальных уравнений:

U T U * = U T y. (4.10) В общем случае, если U T U является невырожденной матрицей, оценки * по методу наименьших квадратов получаются решением матричного уравнения (4.10):

* = [U T U ]1U T y. (4.11) В некоторых случаях, когда матрица U является квадратной мат рицей, что имеет место, например, если размер выборки равен числу оцениваемых параметров, или при использовании регрессионного МНК, и имеет обратную матрицу, то с учетом соотношения [U T U ] 1U T = U 1, вектор оценок может быть определен более про стым способом * = U 1 y. (4.12) Во многих случаях функции достаточно общего вида могут быть разложены в ряд по системе ортонормальных функций:

0, i k T u i (t )u j (t )dt = ij, где ij = - символ Кронекера.

1, i = k В этом случае U T U = I, где I – единичная матрица, и оценки по МНК в базисе ортонормальных функций получаются проще * = U T y.

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

Будем считать, что процедура структурной идентификации вы полнена на предшествующем этапе и порядки числителя и знамена теля передаточной функции модели n и m однозначно заданы. Пусть измерения выполнены на интервале из (n+N) моментов времени и, следовательно, имеются выборки из N измерений для входного и вы u (k ) = [u (0), u (1),..., ( N 1)]T ходного сигналов: и y (k ) = [ y (n), y (n + 1),..., y (n + N )]T. На их основе по каждым k экспе риментально сделанным измерениям входного и выходного сигналов можно приближенно рассчитать (предсказать) следующее k+1 значе ние выходной величины. Такое предсказанное значение можно счи тать его некоторой оценкой, сделанной на основе k предшествующих измерений для последующего k+1 момента времени.

Введем следующие обозначения:

u (k ), y (k ) - экспериментальные данные для входного и выходного воздействий соответственно, полученные в k-ый момент времени;

y (k ) - предсказанное значение выходного сигнала в k момент време ни, рассчитанное по совокупности k-1 предшествующих измерений.

Запишем АРСС-модель идентифицируемого объекта при задан ных порядках n и m. Будем рассматривать объект без запаздывания, т.к. учет запаздывания не вносит принципиальных особенностей в решение задачи и не меняет размерности расширенного вектора дан ных и вектора параметров модели, а лишь приводит к появлению за держки в управляющем сигнале на целое число d периодов квантова ния. Для каждого момента k предсказанное значение выходного сиг нала y (k ) определяется зависимостью:

y (k ) = a1 y (k 1) +... + an y (k n) + b1u (k 1) +... + bmu (k m).

(4.13) На основе (4.13) система соотношений для предсказаний для всей временной выборки из N измерений имеет вид:

u (m ) y (n) y (n 1) u (1) a... y (0)...

u (2)...

y (n + 1)... y (1) u (m + 1)...

y ( n) an......

= b...

y (n + N ) y (n + N 1)... y ( N ) u (m + N ).. u ( N + 1) bm (4.14) Здесь столбец Y T (k ) = [ y (n) y (n + 1)... y (n + N )] представляет вектор предсказанных значений выходного сигнала;

матрица u (m ) u (m 1) y (n 1) y (n 2) u (1)... y (0)...

u ( 2) y (n 1) u (m + 1) y ( n)... y (1) u ( m)...

...

= y (n + N 1) y (n + N 2)... y ( N ) u (m + N ) u (m + N 1)... u ( N + 1) (4.15) представляет определенным образом сформированный массив экспе риментальных данных наблюдений за входным и выходным сигнала T ми;

= [a1... a n b1... bm ] - вектор параметров модели.

В матричной форме соотношение (4.14) имеет вид:

Y (k ) = (k ). (4.16) Разность между векторами измеренных значений выходного сиг нала Y (k ) и предсказанных по модели (4.13) Y (k ) образует ошибку аппроксимации, состоящую из погрешностей измерений выходного сигнала и неточностей значений параметров модели:

e(k, ) = Y (k, ) Y (k, ). (4.17) На основе (4.17) формируется функционал среднеквадратичной ошибки:

n+ N e 2 (k, ).

J ( ) = e ( ) e( ) = T (4.18) k =n Так же, как и ранее, из условия существования минимума J ( ) = 0 определяется выражение для оценки, минимизирую = * щее функцию ошибки:

a1 y ( n)... y (n + 1) an...

= [ T ]1 * T (4.19), b... y (n + N ) bm которое в матричной форме имеет вид:

= [ T ]1 * T Y. (4.20) Полученные выражения (4.19) или (4.20) представляют в явной форме оценку параметров модели методом наименьших квадратов на основе обработки результатов измерений по полной выборке, когда сначала собирается весь объем исходных экспериментальных данных, после чего производится ретроспективная процедура идентификации.

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

Вследствие этого, получили широкое применение рекуррентные вычислительные схемы [19, 39, 57, 58, 71], свободные от указанных недостатков. Сущность рекуррентных процедур состоит в получении оценки вектора параметров (k + 1) на каждом k+1-м шаге путем корректировки оценки на предыдущем k-м шаге. Построение текущей оценки производится на основании k+1 наблюдений и результатов вычисления оценки (k ) на предыдущем шаге схемы [19, 35, 39]:

(k + 1) = (k ) + (k )[ y (k + 1) U (k + 1) (k )], (4.21) где U (k + 1), y (k + 1) - вновь поступающие данные, соответствующие k+1-му наблюдению входного и выходного сигналов;

(k ) - вектор коррекции предыдущей оценки на основании текущих данных, вы числяемый следующим образом:

P(k ) U T (k + 1) (k ) = (4.22), I + U (k + 1) P(k ) U (k + 1) T где I - единичная матрица соответствующей размерности.

Вспомогательный вектор P(k ), содержащий текущие значения входного сигнала, должен быть рассчитан заранее для подготовки к каждому очередному шагу в соответствии с соотношением:

P(k + 1) = [ I (k )U (k + 1)]P(k ). (4.23) Размерность P (k ) не зависит от размера выборки и номера наблюде ния и равна [n n] или [(n + m) (n + m)] при использовании модели (4.5) или (4.16) соответственно.

При использовании АРСС-модели вектор или матрица исходных данных U заменяется на матрицу входо-выходных данных (4.15), и соответствующие формулы рекуррентного МНК приобретают вид:

(k + 1) = (k ) + (k )[ y (k + 1) (k + 1) (k )], (4.24) P(k ) T (k + 1) (k ) = (4.25), I + (k + 1) P(k ) T (k + 1) P(k + 1) = [ I (k ) (k + 1)]P(k ). (4.26) Конструктивная реализация рекуррентного алгоритма вычисле ния вектора параметров на основе МНК сводится к следующим эта пам.

Задаются начальные приближения вектора оценок (0) и вспо могательного вектора P(0). Начальные значения могут быть рассчи таны для некоторого l номера наблюдений на основе стандартной процедуры МНК:

P(l ) = [U lT U l ]1 ;

(l ) = P(l )U lT yl, где U l, yl - выборка из l экспериментальных данных входного и вы ходного сигналов. Можно выбрать (0) произвольно, или используя имеющуюся априорную информацию. Например, можно использо вать следующие значения:

(0) = 0;

P (0) = I, где 1- достаточно большое число.

1. На очередном цикле измерений производится регистрация входного и выходного сигналов и формируется новый вектор данных U T (k + 1), y (k + 1) или (k + 1).

2. Вычисляется вектор коррекции (k ) предыдущей оценки с учетом вновь поступивших данных по формуле (4.22) или (4.25).

3. Определяется вектор новых оценок параметров (k + 1) по формуле (4.21) или (4.24).

4. Производится подготовка к следующему циклу, вычисляется вектор P(k + 1) по формуле (4.23) или (4.26).

Этапы 1 – 4 повторяются на каждом такте процедуры идентифика ции.

Рассмотрим применение полученных алгоритмов реализации МНК для решения некоторых типичных задач идентификации объек тов управления.

4.3 Использование метода наименьших квадратов в задачах идентификации 4.3.1 Идентификация статического объекта регрессионным МНК.

Исходными данными для задачи идентификации являются конеч ный ряд экспериментальных значений входных величин объекта xi и соответствующие значения выходных переменных y j. Математиче ская модель функциональной связи между входными и выходными переменными задается в виде уравнения регрессии [21, 23]:

y M = f (x), (4.27) где f (x) - некоторая аналитическая зависимость, в качестве которой наиболее часто применяются степенные полиномы:

y M = a0 + a1 x + a2 x 2 +... + am x m. (4.28) Задача идентификации ставится как нахождение таких оценок неизвестных параметров ai, i = 0,1,...m, при которых заданная урав нением (4.28) аналитическая зависимость будет наилучшим образом аппроксимировать экспериментальные данные.

В качестве критерия близости используется минимум квадратич ной невязки J значений фактических переменных y j и модельных, рассчитанных по уравнению регрессии (4.28):

( )= N N J (ai ) = e 2 = y j yM j j j =1 j = N = ( y j (a0 + a1 x j + a2 x j +... + am x j )) 2 min, 2 m (4.29) j = где y j – экспериментальное значение выходной переменной, полу ченное в j-ый момент времени;

y M j – модельное (расчётное) значение в тот же момент времени.

Для нахождения коэффициентов регрессии составляют уравнения наличия экстремума по каждому параметру ai :

J = 0, i = 0,1,...m. (4.30) ai Совокупность соотношений (4.30) образует систему уравнений относительно оценок m+1 коэффициентов уравнения регрессии (4.28), решение которой определяет искомые коэффициенты.

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

Пример 4.1 Проиллюстрируем применение метода для решения задачи идентификации в случае аппроксимации опытных данных квадратичным полиномом y M = a0 + a1 x + a 2 x 2 при заданных коэф фициентах регрессии a 0 = 2.5;

a1 = 1.75;

a 2 = 5.06.

Критерий минимума среднеквадратичной ошибки в этом случае определяется функционалом:

[y j (a0 + a1 x j + a 2 x 2j )].

N min J (ai ) = min a0, a1, a2 j = ai Система уравнений для нахождения коэффициентов ai в соответст вии с (4.30) принимает вид:

J ( ) N a = 2 y j a0 a1 x j a2 x j = 0;

j = J ( ) N = 2 y j a0 a1 x j a2 x 2 x j = 0;

(4.31) j a1 j = J ( ) N = 2 y j a0 a1 x j a2 x 2 x 2 = 0.

j j a2 j = Преобразовывая (4.31), получим следующие соотношения:

N N N a0 N + a1 x j + a2 x j = y j ;

j =1 j =1 j = N N N N a0 x j +a1 x j + a2 x j = y j xi ;

2 (4.32) j =1 j =1 j =1 j = N2 N N N a0 x j + a1 x 3 +a2 x 4 = y j x 2.

j j j j =1 j =1 j =1 j = Представим систему (4.32) в матричном виде:

2 N N N xj xj yj N a0 Nj = N j =1 j = N N 3 = y x.

x j j j 1 j j (4.33) x xa j =1 j =1 j =1 j = a2 N N N N x j x j x j y j x j 2 3 4 j =1 j =1 j =1 j = Решением системы (4.33) являются искомые выражения для коэффи циентов уравнения регрессии ai :

N N N xj x2 yj N j a0 N Nj =1 j =1 j = N N a = x 3 y x.

1 j x2 xj j j (4.34) j j =1 j =1 j =1 j = a2 N N N N x2 x 3j xj y j x j 4 j j =1 j =1 j =1 j = В дальнейшем, для удобства использования примем следующие обо значения:

N N N N S1 = x j ;

S 2 = x j ;

S3 = x j ;

S 4 = x j ;

2 3 j =1 j =1 j =1 j = N N N S = y ;

S = y x ;

S = y x 2.

5 j =1 j 6 jj 7 jj j =1 j = В соответствии с принятыми обозначениями, вектор оценок ко эффициентов регрессии ai определяется как решение следующей системы:

a0 n S1 S 2 S a = S S S S. (4.35) 1 1 3 a2 S 2 S 3 S 4 S Приведем программную реализацию рассмотренного метода.

a0=2.5;

% точные коэффициенты регрессии a1=-1.75;

a2=5.06;

N=40;

% размер выборки x=10*normrnd(8, 2, [N 1]);

% моделирование входного воздействия v=0.1*randn(N,1);

% моделирование помехи в виде белого шума y=[a0+a1*x(1:N)+a2*x(1:N).^2+v(1:N)];

% моделирование выходного сиг нала с учетом помехи % формирование по исходным данным суммирующих коэффициентов s1=sum(x(1:N));

s2=sum(x(1:N).^2);

s3=sum(x(1:N).^3);

s4=sum(x(1:N).^4);

s5=sum(y(1:N));

s6=sum(y(1:N).*x(1:N));

s7=sum(y(1:N).*x(1:N).^2);

R=[N s1 s2;

s1 s2 s3;

s2 s3 s4];

%формирование квадратной матрицы дан ных Y=[s5;

s6;

s7];

%формирование вектора данных betta=inv(R)*Y;

% расчет оценок по МНК betta =% рассчитанные оценки параметров 2. -1. 5.0600.

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

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

Положим, что система стационарна и линейна в диапазоне изме нения амплитуды входного сигнала и в окрестностях рабочего режи ма. Исходными данными для идентификации являются эксперимен тальные значения кривой разгона объекта y j, полученные в дискрет ные моменты времени t j, j = 1,2...N.

Применим МНК для определения значений коэффициентов пере даточной функции из условия наилучшего соответствия модели и объекта при установленном заранее на основании формы переходной функции и динамических свойств объекта типе передаточной функ ции.

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

k Wo ( p) =. (4.36) (T1 p + 1)(T2 p + 1) Поставим задачу определения параметров модели объекта k 0, Т 1, Т 2, обеспечивающих наилучшее соответствие модельного описания и экспериментальных данных.

Для решения задачи идентификации перейдем от модели объекта в форме передаточной функции (4.36) к временной характеристике кривой разгона h(t ), вычисляя обратное преобразование Лапласа или с помощью таблиц преобразований h(t ) = L1{h( p )} = L1 Wo ( p ). (4.37) p Для рассматриваемого объекта второго порядка аналитическое представление переходной характеристики имеет вид:

t t k0 T1 T = k 0 1 e T2.

h(t ) = L1 e T1 + (T1 T2 ) p (T1 p + 1)(T2 p + 1) (T1 T2 ) (4.38) Переходную функцию h(t ) будем рассматривать как функцию, параметрически зависящую от коэффициентов k 0, Т 1, Т 2 и от времени t.

Для нахождения параметров переходной характеристики соста вим функционал квадратичной невязки экспериментальных данных yi и расчетных значений по (4.38) для тех же моментов времени t = t j, и минимизируем его:

N t j t j T1 T min J (k, T1, T2 ) = min k0 1 y.

e T1 + e T (T T ) j (T1 T2 ) k,T1,T2 j =1 k,T1,T 1 (4.39) Примечание 1 Если необходимо лишь определить интересующие параметры k 0, Т 1, Т 2 объекта управления, не ставя конкретно задачу выбо ра и обоснования метода решения, скорости сходимости метода, допусти мой точности и т.п., то современные специализированные программные средства для выполнения инженерных и научных расчетов (MathCad, MatLab и подобные) позволяют решить данную задачу, не рассматривая подробно особенности того или иного метода решения, а лишь правильно сформулировав функционал идентификации общего вида (4.6) или, в ча стном случае (4.39), и записав его с учетом синтаксиса среды разработки выбранного программного средства.

Для получения оценок искомых параметров k 0, Т 1, Т 2 методом наименьших квадратов запишем для (4.39) условия минимума функ ционала J k = 0, J = 0, (4.40) T J = 0, T и разрешим полученную систему (4.40). Конкретный способ решения построенной системы приводит к различным вариантам применения МНК.

Решение системы уравнений (4.40), составленной для непрерыв ной модели объекта (4.38) вызывает определенные сложности, свя занные с аналитическими вычислениями. Кроме того, учитывая, что исходными данными являются массивы значений входных и выход ных сигналов, полученные в дискретные моменты времени с опреде ленным периодом квантования t, часто оказывается удобнее при менять дискретные модели объекта вида (4.13).

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

d2y dy T1 2 + T2 + y (t ) = k0u (t ). (4.41) dt dt Процедура дискретизации модели (4.41) приводит к уравнению линейной регрессии, частному случаю общего уравнения (4.13), для которого определены порядки n = 2, m = y (k ) = a1 y (k 1) + a 2 y (k 2) + bu (k 1), k = 1,2,...N, (4.42) где k 0 t 2T1 T2 t (T2 t )t T a1 = ;

a2 = ;

b= (4.43) T1 T1 T - параметры дискретной модели, подлежащие оцениванию;

t период квантования.

4.3.3 Идентификация динамического объекта регрессионным МНК Применим для оценивания параметров модели (4.42) регрессион ную процедуру метода наименьших квадратов.

Пусть накоплено N точек измерения входного и выходного сиг налов объекта. С учетом порядка дискретной модели (n=2), функцио нал, минимизирующий квадратичную ошибку идентификации, будет иметь вид:

N [ y (k ) (a1 y (k 1) + a2 y (k 2) + bu (k 1))] J= min. (4.44) k = Система уравнений для нахождения неизвестных параметров a1, a 2, b, отвечающая равенству нулю соответствующих частных про изводных, имеет вид:

J N = 2 [ y (k ) a1 y (k 1) a2 y (k 2) bu (k 1)]y (k 1) = 0;

a k = J N = 2 [ y (k ) a1 y (k 1) a2 y (k 2) bu (k 1)]y (k 2) = 0;

(4.45) a k = J = 2 N [ y (k ) a y (k 1) a y (k 2) bu (k 1)]u (k 1) = 0.

b 1 k = Проведем ряд преобразований и обозначим суммы произведений соответствующих отдельных измерений в дискретные моменты вре мени t = kt следующим образом:

N N N N S1 = y 2 (k );

S 2 = y (k ) y (k 1);

S3 = y 2 (k 1);

y (k ) y (k 2);

S4 = k =1 k =2 k =2 k = N N N y (k )u (k 1);

S 6 = y ( k 1) y (k 2);

y (k 1)u (k 1);

S5 = S7 = k =2 k =3 k = N N N y y (k 2)u (k 1);

u 2 (k 1).

S8 = (k 2);

S9 = S10 = k =3 k =3 k = Система уравнений (4.45) с учетом принятых обозначений примет вид:

J = 2 S 2 + 2a1S3 + 2a2 S 6 + 2bS 7 = 0;

a J = 2S 4 + 2a1S 6 + 2a2 S8 + 2bS9 = 0;

(4.46) a J = 2 S5 + 2a1S 7 + 2a2 S9 + 2bS9 = 0.

b Представим (4.46) в форме матричного уравнения:

S3 S 6 S 7 a1 S S S S9 a2 = S 4. (4.47) 6 S 7 S9 S10 b S Разрешая (4.47) с помощью регрессионной процедуры МНК (4.12), найдем вектор оцениваемых параметров дискретной модели:

a1 S3 S 6 S 7 S a = S S S9 S 4. (4.48) 2 6 b S 7 S9 S10 S На основе решения системы (4.48) учитывая соотношения (4.43) далее нужно определить параметры непрерывной модели k 0, Т 1, Т 2, однозначным образом связанные с идентифицируемыми параметрами дискретной модели a1, a 2, b :

(a 2 + 1)T1 + t t 2 bT T1 = ;

T2 = ;

k0 = 1.

1 a1 a 2 t t Пример 4.2 Приведем программную реализацию регрессионной процедуры оценивания параметров дискретной и непрерывной моде лей по входным и выходным (незашумленным и зашумленным) дан ным для объекта второго порядка с коэффициентами k 0 = 25;

Т 1 = 36;

Т 2 = 15.

s1=tf([25],[36 15 1])% передаточная функция непрерывной модели объекта T_end=60;

% интервал измерений dt=0.2;

% шаг дискретизации t=0:dt:T_end;

% массив дискретного времени N=length(t);

% размер выборки u=ones(N,1);

% единичное входное воздействие v=0.1*randn(N,1);

%моделирование помехи (при учете) в виде белого шу ма y=lsim(s1,u,t) %+v;

% выходная величина % формирование по исходным данным суммирующих коэффициентов S1=sum(y(1:N).^2);

S2=sum(y(2:N).*y(1:N-1));

S3=sum(y(1:N-1).^2);

S4=sum(y(3:N).*y(1:N-2));

S5=sum(y(2:N).*u(1:N-1));

S6=sum(y(2:N-1).*y(1:N-2));

S7=sum(y(1:N-1).*u(1:N-1));

S8=sum(y(1:N-2).^2);

S9=sum(y(1:N-2).*u(2:N-1));

S10=sum(u(1:N-1).^2);

A=[S3 S6 S7;

S6 S8 S9;

S7 S9 S10];

%формирование квадратной матрицы данных B=[S2 S4 S5]';

%формирование вектора данных betta=inv(A)*B;

% оценки параметров дискретной модели a1= betta(1);

a2= betta(2);

b= betta(3);

T1=dt^2/(1-betta(1)-betta(2));

% расчет параметров непрерывной модели T2=(betta(2)*T1+T1+dt^2)/dt;

K=betta(3)*T1/dt^2;

s2=tf([K],[T1 T2 1]);

y2=lsim(s2,u,t);

plot(t,y,t,y2,':');

% сравнение переходных характеристик объекта и модели grid;

Рассчитанные оценки параметров дискретной модели без учета поме хи:

a1 = 1.9250;

a2 = -0.9260;

b = 0.0247;

Рассчитанные оценки параметров непрерывной модели без учета по мехи T1 = 40.3510;

T2 =15.1222;

K =24.9271;

Экспериментальная переходная характеристика 15 Модельная y(t), ym(t) y(t), ym(t) переходная характеристика 5 0 10 20 30 40 50 0 10 20 30 40 50 Время, с Время, с а б Рисунок 4.


Сравнение переходных характеристик объекта и модели в случае отсутст вия помехи (а) и при учете аддитивной помехи на выходе типа белый шум I=0.1. (б) Сравнение полученных оценок параметров модели с их истинны ми значениями при отсутствии случайных возмущений показывает высокую точность оценивания. Графическое сопоставление (рисунок 4.1) идентифицированных и истинных характеристик объекта пока зывает практически точное совпадение результатов при отсутствии помех и удовлетворительное соответствие выходных сигналов объек та и модели при зашумленных входных данных. Воздействующие на объект помехи существенно влияют на точность оценивания, и их це лесообразно предварительно отфильтровывать. Кроме того, точность полученных оценок зависит от шага дискретизации и выбранного ин тервала измерений.

При использовании МНК получаемые оценки вычисляются с не которыми ошибками, которые называются смещением оценок. Для получения достаточно представительных результатов необходимо выполнить ряд условий [74]:

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

• Объем исследуемой выборки N должен быть достаточным для получения представительных оценок, причем, чем меньше уровень тестового сигнала, тем больше должно быть число N. Целесообразно применять рекуррентный метод наименьших квадратов, который по зволяет в реальном времени получать текущие оценки параметров объекта и по их сходимости определять величину N и момент окон чания эксперимента.

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

4.3.4 Идентификация динамического объекта явным МНК Пример 4.3 Рассмотрим применение явной формы МНК для па раметрической идентификации той же АРСС - модели объекта второ го порядка (4.42), с учетом заданных порядков n = 2, m = 1.

Использование модели (4.42) для оценок коэффициентов a1, a 2, b на основе выборки из N (от 1 до N) экспериментальных данных при водит к следующей системе уравнений вида (4.4):

a1 y (1) + a2 y (0) + bu (1) = y (2);

a y (2) + a y (1) + bu (2) = y (3);

1 (4.49)...

a1 y ( N 1) + a2 y ( N 2) + bu ( N 1) = y ( N ).

Матричная форма записи данной модели имеет стандартный вид (4.5) линейной модели:

y (1) u (1) y ( 2) y ( 0) a y (2) u ( 2) y (3) y (1) a2 =. (4.50)...

...

b y ( N 1) y ( N 2) u ( N 1) y ( N ) С учетом обозначения матрицы исходных входо-выходных данных y (1) u (1) y (0) y (2) u (2) y (1) =, (4.51)... y ( N 1) y ( N 2) u ( N 1) параметры дискретной модели a1, a 2, b определяются на основе об щего соотношения МНК (4.11) следующим образом:

y ( 2) a1 [ ] a = T 1 T y (3). (4.52) 2...

b y ( N ) Приведем программную реализацию явного МНК объекта второ k го порядка с передаточной функцией Wo ( p) = и коэф T1 p 2 + T2 p + фициентами k 0 = 25;

Т 1 = 36;

Т 2 = 15.

s1=tf([25],[36 15 1])% непрерывная передаточная функция объекта T_end=60;

% интервал измерений dt=0.2;

% шаг дискретизации t=0:dt:T_end;

% массив дискретного времени N=length(t);

% размер выборки u=ones(N,1);

% моделирование единичного входного воздействия y=lsim(s1,u,t);

% моделирование выходного воздействия n=2;

% порядок объекта R=[y(n:N-1) y(n-1:N-2) u(n:N-1)];

% формирование расширенной матрицы данных Y=y(n+1:N);

% формирование вектора выходных данных betta=inv(R'*R)*R'*Y;

% расчет параметров непрерывной модели T1=dt^2/(1-betta(1)-betta(2)) T2=(betta(2)*T1+T1+dt^2)/dt K=betta(3)*T1/dt^ Рассчитанные оценки параметров дискретной модели:

a1 = 1.9190;

a2 = -0.9200;

b = 0.0266;

Рассчитанные оценки параметров непрерывной модели T1 = 37.5243;

T2 =15.2014;

K = 25.0000;

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

4.3.5 Идентификация динамического объекта рекуррентным МНК Пример 4.4 Приведем программную реализацию оценивания параметров k 0, Т 1, Т 2 объекта из предыдущего примера с помощью рекуррентного МНК (4.24) при использовании АРСС - модели объек та второго порядка (4.42).

s1=tf([25],[36 15 1]) % непрерывная передаточная функция объекта T_end=60;

% интервал измерений dt=0.2;

% шаг дискретизации t=0:dt:T_end;

% массив дискретного времени N=length(t);

% размер выборки u=ones(N,1);

% массив значений единичного входного воздействия y=lsim(s1,u,t);

%массив значений выходного воздействия n=2;

% порядок объекта I=diag([1 1 1]);

i=1;

% начальный шаг P=1000*I;

% начальное приближение betta=[0;

0;

0];

bet(i,:)=betta;

% массив оценок параметров % очередной шаг вычислений for i=n:N- R=[y(i+n-2:-1:i-1);

u(i+n-2:-1:i)]';

% формирование расширенной матрицы данных gamma=P*R'/(R*P*R'+1);

betta=betta+gamma*(y(i+1)-R*betta);

P=(I-gamma*R)*P;

bet(i,:)=betta;

end;

plot(bet,'+');

T1=dt^2/(1-betta(1)-betta(2)) % расчет параметров непрерывной модели T2=(betta(2)*T1+T1+dt^2)/dt K=betta(3)*T1/dt^ Оценки параметров непрерывной модели:

T1 = 35.6366;

T2 =15.4333;

K = 25.0975.

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

bet 1. bet 1(k), bet 2(k), bet 3(k) 0. bet -0. bet - 0 50 100 150 200 250 k Рисунок 4. Сходимость оценок параметров дискретной модели 4.3.6 Определение импульсной переходной функции объекта с помощью метода наименьших квадратов Рассмотрим использование МНК для идентификации импульсной переходной функции (ИПФ) линейного стационарного объекта с одним входом и одним выходом. В соответствии с рассмотренной ранее схемой проведения эксперимента (рисунок 3.7), требуется определить ИПФ по результатам измерений входного u (t ) и выходного сигналов на конечном промежутке времени y (t ) длительностью Т в условиях действия помехи (t ) типа белого шума, приведенной к выходу.

Выходной сигнал линейной стационарной системы при нулевых начальных условиях выражается стандартным интегралом свертки:

T y (t ) = w(t )u (t )d + (t ), (4.53) где w(t ) – импульсная переходная функция.

Проведем временную дискретизацию уравнения (4.53), с равно мерным интервалом квантования t. Выходной сигнал в произволь ный момент времени t = jt определяется следующим соотношени ем:

N s w(i)u[( j i)t ]t + j, y ( j ) = j = 0,1,... N m 1, (4.54) i = где Tm = N m – время измерения выходного сигнала;

Ts = N s - вре мя оценивания, т.е. установления реакции ИПФ (не более 5% от сво его пикового значения).

Запишем выражение (4.54) в компактном виде:

N s wi u j i t + j, yj = j = 0,1,... N m 1, (4.55) i = где wi = w(i ) ;

u j i = u[( j i )t ];

y j = y ( j ).

Величина j содержит как невязку в дискретные моменты време ни ( j ), так и ошибки, возникающие за счет аппроксимации непре u (t ) рывной зависимости кусочно-постоянной функцией u[( j i )].

Проведенная процедура дискретизации во времени (4.54) приво дит к тому, что оценивание непрерывной функции w(t ) заменяется оцениванием конечного множества параметров w0,...wN s 1.

Выражения (4.55) в развернутом виде представляются следую щим образом:

... u ( N s 1) w0 t y0 u 0 u y u... u ( N s 2 ) w1t u = +............

u N m N s wN s 1t N m y N m 1 u N m 1 u N m (4.56) или в матричной форме:

y = U +, (4.57) где = w t - вектор-столбец идентифицируемых параметров, y, и U –вектор-столбцы и матрица соответствующих выборочных значений.

Таким образом, оценивание ИПФ сводится к оцениванию вектора параметров при заданной матрице U и векторе измерений у. Ре зультатом оценивания является нахождение вектора, минимизи рующего сумму квадратов невязок на интервале измерения:

J ( ) = ( y U )T ( y U ) min. (4.58) Оценка по МНК * удовлетворяет требованиям J = min J = J = * и находится из условия экстремума функционала (4.58):

J = 2U T ( y U * ) = 0. (4.59) * = Система уравнений (4.59) в матричной форме имеет вид:

U T U * = U T y, (4.60) и ее решение относительно вектора параметров находится следую щим образом:

[ ] * = U T U U T y. (4.61) Соответственно, выражение явной формы метода наименьших квад ратов (4.11) для оценивания конечного множества параметров им пульсной переходной характеристики принимает следующий вид:

[ ] 1 T 1 T w= (4.62) U U U y.

t Перепишем уравнение (4.60) относительно сумм выборочных значе ний сигналов:

N m 1 N s 1 N m u j i u j k wi* t = u j i y j, Nm N m j = i =0 j = i = 0,1,... N s 1;

k = 0,1,... N s 1 (4.63) Уравнение (4.63) в непрерывной форме представляет известное урав нение Винера-Хопфа 1 m * Ts T T 1m T u (t )u (t )dt w ( )d = T u (t ) y (t )dt, (4.64) 0 m 0 m Tm K uu ( )w ( )d = K uy ( ), * или (4.65) 1 Tm где K uu ( ) = u (t )u (t )dt - автокорреляционная функция вход Tm 1 Tm ного сигнала;

K uy ( ) = u (t ) y (t )dt - взаимная корреляционная Tm функция входного и выходного сигнала.

Пример 4.5 Приведем программную реализацию процедуры идентификации импульсной переходной характеристики для того же k базового объекта с передаточной функцией Wo ( p) = и T1 p + T2 p + коэффициентами k 0 = 25;


Т 1 = 36;

Т 2 = 15.

s1=tf([25],[36 15 1])% непрерывная передаточная функция объекта T_end=45;

% интервал измерений dt=1.5;

% шаг дискретизации t=0:dt:T_end;

% массив дискретного времени N=length(t);

% размер выборки u=sign(normrnd(0, 2, [N 1]));

%моделирование входного воздействия y=lsim(s1,u,t) ;

% выходное воздействие % заполнение матрицы U входных значений for i=1:N for j=1:N if(i=j) U(i,j)=u(i+1-j);

else U(i,j)=0;

end;

end end w=1/dt*(inv(U'*U))*(U'*y);

w0= impulse(s1,t);

plot (t, w0, t, w,':');

grid;

Результаты идентификации ИПФ рассматриваемого объекта представлены на рисунке 4.3.

1. 1. 0. w(t), w0(t) 0. 0. 0. 0 5 10 15 20 25 30 35 40 Время, с Рисунок 4. Аналитическая (1) и экспериментальная (2) импульсные характеристики объекта Идентификационный эксперимент проводился при подаче на объект случайной последовательности сигналов ± 1. Сравнение ана литической и экспериментально полученной импульсных весовых характеристик (рисунок 4.3) показывает высокую точность данного метода, существенно зависящую от интервала измерений и шага дис кретизации. Расчеты показывают, что нужно внимательно выбирать интервал измерений Tm и время установления реакции Ts = N s t, т.к.

при проведении процедуры оценивания на участке, где ИПФ стре [ ] мится к нулю, матрица U T U становится близкой к вырожденной, что влечет расходимость решения.

4.4 Градиентные методы Рассмотрим общую задачу минимизации квадратичной невязки выходных сигналов модели и объекта для функционала t J = [ yo (t ) yM (t )] dt.

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

(k + 1) = (k ) + (k ) gradJ [ (k )], (4.67) где (k ) - текущее приближение к истинному вектору параметров * ;

(k ) - служебный параметр, характеризующий длину k-го шага итерационного процесса;

k - номер итерации.

Для определения направления движения к экстремуму использу ется градиент – n- мерный вектор, составляющие которого являются частными производными функции f (x), вычисленными в точке х:

J J J J ( ) = (4.68),,....

1 2 n Градиент указывает направление наискорейшего роста функции (об ратное направление будет направлением наискорейшего спуска).

Градиент функции может быть определен аналитически, а если функ ция J ( ) не задана, то с помощью экспериментов.

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

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

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

Траектория движения в каждой точке ортогональна к линиям уровня J ( ) = const. Итерационная процедура определения значения вектора параметров на очередном шаге имеет вид:

(k + 1) = (k ) + (k ) s (k ), (4.69) где s (k ) - вектор, определяющий направление спуска, находится сле дующим образом:

J ( (k )) s(k ) =. (4.70) J ( (k )) Нормирующий коэффициент вектора градиента целевой функции определяется соотношением:

J [ n (k )] J [1 (k )] J [ 2 (k )] 2 J [ (k )] = + +... +.

+ (4.71) x x x 1 2 n Величина шага (k ) выбирается таким образом, чтобы целевая функция в выбранном направлении не перестала убывать. Значение (k ), в соответствии с этим, находится из условия минимума квадра тической аппроксимации целевой функции J ( ) по (k ) в точке (k ) :

[J ( (k ))]T s(k ), (k ) = (4.72) [s(k )]T 2 J ( (k )) s(k ) где 2 J ( (k )) - матрица вторых производных:

2J 2 J 2J...

1 2 1 n 2J 2 J 2J...

2 J ( (k )) = (4.73) 2 n.

2............

2J 2 J 2J...

n 1 n 2 n В точке (k + 1) определяется новое направление движения к минимуму, которое будет ортогонально предыдущему и далее повто ряется процесс приближения к экстремуму.

Для реализации алгоритма данного метода нужно задать точку начального приближения (0) и последовательно на каждом шаге алгоритма k = 0,1,... вычислять следующие компоненты:

1) вектор градиента целевой функции J ( (k )) в точке = (k ) соответственно (4.68);

нормирующий коэффициент вектора градиента J [ (k )] 2) соответственно (4.71);

3) вектор s (k ) соответственно (4.70);

матрица вторых производных 2 J ( (k )) соответственно 4) (4.73);

значение шага (k ) соответственно (4.72);

5) (k + 1) соответственно 6) новое значение приближения (4.69).

Градиентные методы являются основой для идентификации дос таточно сложных объектов, для оптимизации нелинейных критериев качества идентификации.

Пример 4.6 Рассмотрим задачу идентификации инерционного объекта второго порядка, описываемого передаточной функцией ko Wo ( p) =. Требуется найти параметры модели (T1 p + 1) (T2 p + 1) km пониженного порядка, имеющей вид Wm ( p) =, такие, чтобы Tm p + квадратичный критерий t t t t кон T1 T k dt (4.74) k 1 e 1 Tm T1 T F= e + e o m ( T2 + T1 ) ( T2 + T1 ) принимал свое минимальное значение при известных параметрах объекта k0 = 27.5;

T1 = 11.8;

T2 = 3.2. Примем период времени, на кото ром проводится идентификация параметров t кон = 100с.

В такой постановке задачи вектор искомых параметров объекта k m содержит два параметра =. После интегрирования и выполне Tm ния необходимых преобразований в (4.74) получают явную форму функционала качества:

2 3k mTm k F (k m, Tm ) = (ko k m ) t кон + 2(T1 T2 ) 2 (T1 + T2 ) 2k 0 k m (3T14 3T24 + 4T12T22 + T1T2 (T12 + T22 )) + (T1 + Tm )(T2 + Tm )(T1 T2 ) (T1T2 (T12 T22 ) + (T1T2 (T1 T2 ) + T13 T23 )Tm + (T12 T22 )Tm + (T1 T2 )Tm ) 2 (4.75) Обозначим коэффициент, не зависящий от искомых параметров, k как с = (3T14 3T24 + 4T12T22 + T1T2 (T12 + T22 )).

2(T1 T2 ) (T1 + T2 ) Коэффициенты при различных степенях искомого параметра Tm обо значим следующим образом:

d 0 = T1T2 (T12 T22 );

d1 = T1T2 (T1 T2 ) + T13 T23 ;

d 2 = T12 T22 ;

d 3 = T1 T2.

С учетом принятых обозначений выражение (4.75) запишем следую щим образом:

3k mTm 2k 0 k m F (k m, Tm ) = (k o k m ) t кон + (T1 + Tm )(T2 + Tm )(T1 T2 ) (4.76) 2 (d 0 + d1Tm + d 2Tm + d 3Tm ) + с.

На основе (4.76) в точках (k ), k = 1,2,... вычисляются компоненты вектора градиента - первые производные критерия идентификации по искомым параметрам:

F 2k 0 2 = 2 t кон ( k o k m ) 3k mTm + (d 0 + d1Tm + d 2Tm + d 3Tm );

k m (T1 + Tm )(T2 + Tm )(T1 T2 ) 3k F 2k 0 k m = m + ((T1 + Tm )(T2 + Tm ) Tm (T1 + Tm ) 2 (T2 + Tm ) 2 (T1 T2 ) );

(d1 + 2d 2Tm + 3d 3Tm ) (T1 + T2 + 2Tm )(d 0 + d1Tm + d 2Tm + d 3Tm ) 2 2 T F F и составляется вектор градиента F (k m, Tm ) =.

Tm k m Далее определяется нормирующий коэффициент 2 F F F (k m, Tm ) = k + T. Затем составляется матрица вторых m m 2F 2F k m Tm k m производных F (k m, Tm ) =, компоненты кото F 2F T k Tm mm рой определяются следующим образом:

2F = 2 t кон 3Tm ;

k m 2F 2k 0 k m = ((T1 + Tm ) 2 (T2 + Tm ) 2 ((T1 + Tm )(T2 + Tm ) (T1 + Tm ) (T2 + Tm ) (T1 T2 ) 4 Tm ) (2d 2 + 6d 3Tm ) 2(d 0 + d1Tm + d 2Tm + d 3Tm ) 2(T1 + Tm )(T2 + Tm ) 2 (T1 + T2 + 2Tm )((T1 + Tm )(T2 + Tm )(d1 + 2d 2Tm + 3d 3Tm ) )) (T1 + T2 + 2Tm )(d 0 + d1Tm + d 2Tm + d 3Tm ) 2 2F 2k 0 k m = ((T1 + Tm ) 2 (T2 + Tm ) 2 ((T1 + Tm )(T2 + Tm ) 2 4 (T1 + Tm ) (T2 + Tm ) (T1 T2 ) Tm ) (2d 2 + 6d 3Tm ) 2(d 0 + d1Tm + d 2Tm + d 3Tm ) 2(T1 + Tm )(T2 + Tm ) 2 (T1 + T2 + 2Tm )((T1 + Tm )(T2 + Tm )(d1 + 2d 2Tm + 3d 3Tm ) )) (T1 + T2 + 2Tm )(d 0 + d1Tm + d 2Tm + d 3Tm ) 2 2F 2F 2k = = 3k m + ((T1 + Tm )(T2 + Tm ) Tm k m k m Tm (T1 + Tm ) (T2 + Tm ) 2 (T1 T2 ) ) 2 2 (d1 + 2d 2Tm + 3d 3Tm ) (T1 + T2 + 2Tm )(d 0 + d1Tm + d 2Tm + d 3Tm ) Вычислительная процедура заканчивается по достижении норми рованной разности двух соседних приближений некоторого заданно го числа eps=0.50.

(k + 1) (k ) eps.

(k ) Рассмотрим пример программной реализации итерационной проце дуры.

k0=27.5;

% Задание параметров объекта T1=11.8;

T2=3.2;

km=1;

% Задание начальных приближений искомых параметров Tm=1;

z(:,1)=[ 1;

1];

z(:,2)=[ 1;

1];

i=2;

while (abs(z(i)-z(i-1))/abs(z(i))0.5) || (i5) i=i+1;

x=[km;

Tm];

% вектор параметров c0=k0^2*(-3*T1^4 3*T2^4+4*T1^2*T2^2+T1*T2*(T1^2+T2^2))/(2*(T1 T2)^2*(T1+T2));

c1=(k0-km)^2*100;

c2=-3*km^2*Tm/2;

b1=2*k0 /((T1+Tm)*(T2+Tm)*(T1-T2));

b2=(Tm*(T1-T2)+T1^2 T2^2)*(T1+Tm)*(T2+Tm)+Tm*(T2^2*(T1+Tm)-T1^2*(T2+Tm));

c3=b1*b2;

d0= T1*T2*(T1^2-T2^2);

d1= (T1*T2*(T1-T2)+ T1^3-T2^3);

d2=(T1^2-T2^2) ;

d3=(T1-T2) ;

D0= d0+d1*Tm+d2*Tm^2+d3*Tm^ F=c1+c2+b1*km *D0+c0;

% функционал качества D1=d1+2*d2*Tm+3*d3*Tm^ % первые производные функционала f1=-2*(k0-km)*100-3*km*Tm+b1*D0;

f2=-3*km^2/2+b1*km/((T1+Tm)*(T2+Tm))*(D1*(T1+Tm)* (T2+Tm)-D0* (T1+T2+2*Tm));

gr=[f1;

f2];

% вектор градиента nor_gr=(f1^2+f2^2)^(1/2);

% нормирующий коэффициент s=-gr/nor_gr;

% вектор s % вторые производные функционала df_km_2=2*100-3*Tm;

R=(T1+Tm)^2*(T2+Tm)^2*((T1+Tm)*(T2+Tm)*(2*d2+6*d3*Tm) 2*D0)-2*(T1+Tm)*(T2+Tm)*(T1+T2+2*Tm)* ((T1+Tm)*(T2+Tm)*D1-(T1+T2+2*Tm)* D0);

df_Tm_2= b1*km /((T1+Tm)^3*(T2+Tm)^3)*R;

df2_km_Tm=-3*km+b1/((T1+Tm)*(T2+Tm))*(D1*(T1+Tm)* (T2+Tm)- D0*(T1+T2+2*Tm));

gr2=[df_km_2 df2_km_Tm;

df2_km_Tm df_Tm_2];

% матрица вто рых производных l=-gr'*s/(s'*gr2*s);

% величина шага x=x+l*s;

% новое значение вектора параметров km=x(1);

Tm=x(2);

z(:,i)=x;

end;

figure(1);

plot(z(1,3:i), z(2,3:i));

grid;

s1=tf([k0],[T1*T2 (T1+T2) 1]) % передаточная функция объекта s2=tf([km],[Tm 1]) % передаточная функция модели t=0:0.1:100;

N=length(t);

u=ones(N,1);

y1=lsim(s1,u,t);

% переходная функция объекта y2=lsim(s2,u,t);

% переходная функция модели figure(2);

plot (t, y1, t, y2, ':b');

grid;

Реализация итерационной процедуры до достижения заданной точности дает следующие параметры модели:

k m = 27.8805;

Tm = 15.5663. Значение функционала качества: F = 64.8810.

На рисунках 4.4, 4.5 приведены графические результаты расчета.

Tm 23.5 24 24.5 25 25.5 26 26.5 27 27.5 km Рисунок 4. Сходимость параметров модели в двухпараметрическом пространстве Переходная характеристика модели y(t), y m(t) Переходная характеристика объекта 0 10 20 30 40 50 60 70 80 90 Время, с Рисунок 4. Переходные характеристики объекта и модели На рисунках 4.4 и 4.5 приведены траектория сходимости итера ционной процедуры поиска параметров k m, Tm и результаты сопос тавления истинной и идентифицированной переходных характери стик, показывающие удовлетворительное качество идентификации.

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

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

Кроме того, реализация градиентных алгоритмов требует определен ного объема аналитических вычислений, что может вызывать значи тельные сложности. Рекомендуется совмещать применение различ ных итерационных методов с учетом свойств каждого из них [74].

5 Оценивание состояния объекта Зная прошлое, легко предсказать будущее, но изменить его невозможно Из концепции детерминизма 5.1 Общий подход к задаче оценивания переменных состояния В настоящее время широкое распространение получили методы оценивания (идентификации) в пространстве состояний, основанные на работах Р. Калмана и Л. Заде. Эти методы позволяют в реальном времени получать алгоритмы обработки данных с использованием современного уровня развития средств вычислительной техники. Раз личные варианты алгоритмов обработки результатов измерений при меняются при решении задач оценивания состояний в системах управления, оптимальных, самонастраивающихся и адаптивных сис темах, измерительных комплексах и средствах, теории случайных сигналов, навигационной аппаратуре, медицине и многих других от раслях.

На практике достаточно распространенной является ситуация, ко гда не все компоненты вектора состояний доступны для измерения. В этом случае, чтобы в системе управления возможно было использо вать обратную связь по состоянию, необходимо восстановить вектор состояния системы, недоступный для измерения. Восстановление вектора состояния называется его оценкой, а устройства, формирую щие на выходе вектор оценки состояний, а также позволяющие отде лить полезный сигнал от помех, – наблюдателями (идентификатора ми, фильтрами) [5, 6, 19, 31, 35, 39, 51, 57]. Алгоритмы работы на блюдателей обычно строятся таким образом, чтобы при обработке сигналов получалась оптимальная в определенном смысле оценка не которой переменной. Фильтры применяются для следующих проце дур:

• интерполяции, сглаживания или экстраполяции сигналов;

• оценивания состояния объектов;

• оценивания параметров объектов и т.д.

Рассмотрим стандартную ситуацию, при которой возникает зада ча оценивания состояний. Положим, что имеется некоторая система, состояние которой в любой момент времени однозначно характеризу ется определённым набором координат состояний (например, про странственными координатами, скоростью, уровнями напряжения и т.д.), т.е. эти величины являются элементами вектора состояния сис темы x(t ) в каждый момент времени t. Часть этих координат, как правило, недоступна для непосредственного определения. Положим далее, что имеется ряд переменных, некоторым образом связанных с состоянием системы, и их можно измерить с заданной точностью, т.е.

эти величины составляют вектор измерений y (t ), относящийся к оп ределённому моменту времени t. Необходимо оценить значения век тора состояний x(t ), недоступного для непосредственного измерения.

5.2 Оптимальный наблюдатель полного порядка (фильтр Калма на) Одной из наиболее употребительных схем наблюдателей в конту рах автоматического управления и регулирования является адаптив ный фильтр рекурсивного типа, известный как фильтр Калмана (Кал мана-Бьюси) [5, 14, 19, 30, 35, 39, 48, 51, 57, 58, 59, 61, 68, 69, 75].

Алгоритм фильтра Калмана позволяет в реальном времени по ^ строить оптимальную оценку состояния системы x(t ), основываясь на измерениях y (t ), неизбежно содержащих погрешности. Вектор измерений y (t ) рассматривается в качестве многомерного выходного сигнала системы, при этом зашумлённого, а вектор состояния x(t ) неизвестный многомерный сигнал, подлежащий определению. Усло вием оптимальности построенной оценки является минимум средне квадратичного отклонения оцененного значения от истинного.

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

1) Определяется n - мерное пространство функций (в общем случае – гильбертово) и далее полагается, что все временные процес сы принадлежат этому пространству.

2) Имеется линейная динамическая стационарная система (фор мирующий фильтр), описываемая уравнением:

dx = Ax (t ) + Bu (t ) + v(t ), (5.1) dt где, кроме ранее рассмотренных обозначений, принятых при по становке задачи в пространстве состояний (1.20), используются сле дующие:

x(t ) – случайный марковский n - мерный процесс, задаваемый необходимой априорной информацией:

M {x} = x0 ;

cov{x} = P0 ;

(5.2) u (t ) – измеряемое векторное входное воздействие, которое может быть как детерминированной, так и случайной величиной;

v(t ) – k -мерный вектор случайных воздействий, полагаемых про цессами типа белого шума:

M {v} = 0;

{ } cov{v} = M v(t )v( )T = V (t ), (5.3) где (t ) - дельта-функция;

V - симметричная неотрицательно опре деленная матрица интенсивности белого шума v(t ) размерностью [k k ] ;

3) Процесс x(t ) наблюдается с помощью измерителя, и вектор измеряемых координат определяется соотношением:

y (t ) = Cx (t ) + (t ), (5.4) где (t ) – p - мерный вектор шумов измерения, полагаемый случай ным процессом в виде белого шума:

M { } = 0;

{ } cov{ } = M (t ) ( )T = R (t ), (5.5) где R - симметричная неотрицательно определенная матрица интен сивности белого шума (t ) размерностью [ p p] ;

4) Процессы v(t ) и (t ), а также x(t ) и v(t ), x(t ) и (t ) полага ются некоррелированными:

{ } { } { } M v(t ) ( )T = 0;

M x(t )v ( )T = 0;

M x(t ) ( )T = 0, формирующий фильтр удовлетворяет условиям физической реали зуемости w(t ) = 0, t 0, а динамическая система, соответствующая (5.1), (5.4), является полностью управляемой и полностью наблюдае мой.

5) Требуется построить линейную динамическую систему, обес печивающую получение оптимальной оценки x(t ) вектора x(t ), если ошибка оценивания задана:

e(t ) = x(t ) x(t ), (5.6) и критерием оптимальности является условие минимума ее квадра тичной нормы:

J = e 2 (t )dt = e(t ) min. (5.7) Задача фильтрации по Калману-Бьюси соответствует следующей структурной схеме:

Рисунок 5. Структурная схема формирующего фильтра и измерителя при действии случайных возмущений и наличии помех ФФ – формирующий фильтр При построении фильтра Калмана используется идея n-мерного наблюдателя (наблюдателя полного порядка), когда в качестве иден тификатора состояния принимается математическая модель системы.

Фильтр Калмана осуществляет процедуру рекурсивного оценива ния на основе наблюдений за входным и выходным сигналами объек та, где для уменьшения дисперсии оценок в алгоритм идентификато ра вводится корректирующая обратная связь по выходу системы y (t ).

Исходя из требования получения несмещенной оценки x(t ), уравнение фильтра имеет вид [5, 74]:

dx ^ = Ax(t ) + Bu (t ) + L y (t ) С x(t );

(5.8) dt x (t 0 ) = x 0, где L - матрица коэффициентов усиления фильтра размером [n 1], обеспечивающая оптимальную в смысле минимальной дисперсии оценку состояния, которая определяется выражением L = PC T R 1. (5.9) В (5.9) P - ковариационная матрица ошибок оценивания, которая в случае стационарности процессов определяется решением алгеб раического матричного уравнения Риккати:



Pages:     | 1 || 3 |
 





 
© 2013 www.libed.ru - «Бесплатная библиотека научно-практических конференций»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.