ОТУ_Поляков_Часть2_2009. Теория автоматического управления для чайников Часть ii. Управление при случайных
 Скачать 1.41 Mb.
|
|
свойств: 1) ) 0 ( X R – это средний квадрат случайного процесса, поэтому всегда 0 ) 0 ( ≥ X R ; для цен- трированных процессов (с нулевым средним) эта величина совпадает с дисперсией; 2) при 0 = τ корреляционная функция имеет наибольшее значение, в том числе и наиболь- шее по модулю, то есть ) 0 ( ) ( X X R R ≤ τ при всех τ ; 3) ) ( ) ( τ τ − = X X R R , то есть ) ( τ X R – симметричная (четная) функция, это доказывается под- становкой τ − вместо τ в интеграл (1); поэтому можно считать корреляционную функ- цию только для 0 ≥ τ , а вторую часть строить симметрично. В качестве примера приведем корреляционную функцию дискретного сигнала, который пере- ключается между значениями A и A − в случайные моменты времени: 0 t ) (t x A A − 0 τ ) ( τ X R 2 A τ α τ − = e A R X 2 ) ( © К.Ю. Поляков, 2009 14 Корреляционная функция имеет вид τ α τ − = e A R X 2 ) ( , где α – среднее число переключений за 1 с. Заметим, что одна и та же корреляционная функция может соответствовать многих совер- шенно различным процессам. Корреляционная функция не всегда положительна. На следующих рисунках показано из- менение ординаты поверхности морского волнения и корреляционная функция этого сигнала (одна из теоретических моделей): Здесь r D – дисперсия волновой ординаты, α – коэффициент затухания и β – средняя частота волнения. Заметим также, что чаще всего корреляционная функция убывает по модулю, т.е. чем дальше от нуля, тем меньше значение модуля корреляционной функции (чем больше расстоя- ние между отсчетами, тем меньше связь между ними). Это справедливо не для всех случайных процессов, но для большинства практических ситуаций. 2.5. Спектральная плотность В теории управления существуют и взаимно дополняют друг друга два подхода: 1) временнóй – исследование процессов во времени; 2) частотный – исследование частотных свойств сигналов и систем (с помощью передаточ- ных функций и частотных характеристик). Аналогичная ситуация наблюдается и при рассмотрении случайных процессов. Основная временная характеристика стационарного процесса – это корреляционная функция, а частотные свойства описываются спектральной плотностью. Спектральная плотность – это функция, которая показывает распределение мощности сигнала по частотам. Такая информация о полезных сигналах, помехах и возмущениях очень важна для разработчика систем управления. Система должна быть спроектирована так, чтобы усиливать сигналы с «полезными» частотами и подавлять «вредные» частоты, характерные для помех и возмущений. Для перехода от временнóго описания детерминированных (не случайных) процессов к частотному, используют преобразования Фурье и Лапласа. Аналогично спектральная плотность случайного процесса может быть найдена как преобразование Фурье от корреляционной функ- ции 4 : { } ) ( ) ( ) ( τ τ τ ω ωτ X j X X R d e R S F = = ∫ ∞ ∞ − − 4 Эта формула называется формулой Винера-Хинчина. Строго говоря, это не определение спектральной плотности, а следствие из него. Математически корректное определение можно найти в литературе [1,3]. 0 t ) (t x 0 τ ) ( τ X R βτ τ τ α cos ) ( − = e D R r X © К.Ю. Поляков, 2009 15 Здесь 1 − = j – мнимая единица, а ω – угловая частота в рад/с ( f π ω 2 = , где f – «обычная» частота в герцах). Используя формулу Эйлера, можно представить экспоненту в виде сумму вещественной (косинусной) и мнимой (синусной) составляющих: ωτ ωτ ωτ sin cos j e j − = − Функция ωτ τ sin ) ( X R – нечетная по τ , поэтому интеграл от нее в симметричных пределах ра- вен нулю. Напротив, функция ωτ τ cos ) ( X R – четная, так что при интегрировании можно взять интервал от 0 до ∞ и удвоить результат: ∫ ∞ = 0 cos ) ( 2 ) ( τ ωτ τ ω d R S X X . (2) Спектральная плотность чем-то похожа на плотность распределения вероятностей, только она характеризует плотность распределения мощности сигнала по частотам. Если случайный процесс – это напряжение в вольтах, то его корреляционная функция измеряется в В 2 , а спек- тральная плотность – в В 2 /Гц. Спектральная плотность случайного процесса, имеющего корреляционную функцию τ α τ − = e A R X 2 ) ( , вычисляется как ∫ ∫ ∫ ∞ + − ∞ − − ∞ ∞ − − − + = = 0 ) ( 2 0 ) ( 2 2 ) ( τ τ τ ω τ ω α τ ω α ωτ τ α d e A d e A d e e A S j j j X Интервал интегрирования разбит на две части. При 0 < τ имеем τ τ − = , а при 0 > τ – τ τ = . Выполняя интегрирование, получаем 2 2 2 2 0 ) ( 2 0 ) ( 2 2 1 1 ) ( α ω α ω α ω α ω α ω α ω τ ω α τ ω α + = ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ + + − = + − − = ∞ + − ∞ − − A j j A j e A j e A S j j X На рисунке слева показана корреляционная функция, а справа – соответствующая ей спек- тральная плотность мощности: Свойства спектральной плотности: 1) это неотрицательная, четная функция угловой частоты ω (график расположен выше оси абсцисс и симметричен относительно вертикальной оси); 2) интеграл от ) ( ω X S на некотором интервале частот ] ; [ 2 1 ω ω дает мощность, которая связана с этими частотами; поскольку функция ) ( ω X S – четная, результат интегрирования на ] ; [ 2 1 ω ω нужно удвоить, чтобы учесть также и полосу ] ; [ 1 2 ω ω − − ; 3) площадь под кривой определяет средний квадрат случайного процесса (для центрирован- ного процесса он равен дисперсии): ∫ ∞ ∞ − = ω ω π d S x X ) ( 2 1 2 0 ω ) ( ω X S 2 2 2 2 ) ( α ω α ω + = A S X 0 τ ) ( τ X R 2 A τ α τ − = e A R X 2 ) ( © К.Ю. Поляков, 2009 16 Множитель π 2 / 1 нужен для согласования единиц измерения, поскольку угловая частота f π ω 2 = измеряется не в герцах, а в рад/с. Учитывая, что функция ) ( ω X S четная, можно интег- рировать ее только при 0 > ω , а результат удвоить: ∫ ∞ = 0 2 ) ( 1 ω ω π d S x X В теории управления нередко записывают спектральную плотность как функцию ком- плексной переменной s , связанной с угловой частотой по формуле ω j s = (отсюда следует 2 2 ω − = s ). Хотя это не совсем корректно с точки зрения математики, мы будем использовать запись ) (s S X для обозначения спектральной плотности ) ( ω X S , в которой выполнена замена 2 2 ω − = s : 2 2 2 2 ) ( α ω α ω + = A S X ⇔ 2 2 2 2 ) ( α α + − = s A s S X 2.6. Гармонический сигнал Рассмотрим гармонический сигнал ( ) θ ω + = t A t x 0 sin ) ( , где θ – случайная фаза, равномерно распределенная в интервале от 0 до π 2 . Три реализации этого процесса (с разными фазами θ ) показаны на рисунке: Это тоже случайный процесс, однако его отличие от «классических» случайных процессов со- стоит в том, что зная (или определив) случайную фазу θ , мы может вычислить значение этого сигнала при любом t . Таким процессы называют квазидетерминированными. Как только фаза θ определена, процесс становится детерминированным (не случайным). Использование формулы для усреднения по времени (1) дает dt t t T A R T T T X ) sin( ) sin( 2 1 lim ) ( 0 0 0 2 τ ω θ ω θ ω τ + + + = ∫ − →∞ После несложных преобразований (применение формулы произведения синусов и интегриро- вание), получим τ ω τ 0 2 cos 2 ) ( A R X = Чтобы найти спектр такого сигнала, вычислим преобразование Фурье для корреляционной функции. По таблицам находим { } [ ] ) ( ) ( cos 0 0 0 ω ω δ ω ω δ π τ ω − + + = F , поэтому [ ] ) ( ) ( 2 ) ( 0 0 2 ω ω δ ω ω δ π ω − + + = A S X 0 t ) (t x © К.Ю. Поляков, 2009 17 Это значит, что спектральная плотность состоит из двух дельта-функций для частот 0 ω − и 0 ω , а в остальных точках равна нулю. Действительно, с самого начала было легко догадаться, что вся энергия такого сигнала сосредоточена на одной частоте. 2.7. Белый шум В математике для теоретических исследований иногда удобно использовать математиче- ские объекты, которые нереализуемы на практике (например, дельта-функцию). В теории слу- чайных процессов важную роль играет белый шум 5 , имеющий равномерную спектральную плотность по всем частотам, то есть, const ) ( 0 = = S S X ω . Очевидно, что при этом площадь под кривой спектральной плотности (определяющая средний квадрат процесса) бесконечна, то есть сигнал имеет бесконечную мощность и не может существовать в природе. Если нет никакой информации о свойствах случайных возмущений, действующих на сис- темы, часто считают, что они приближенно описываются моделью белого шума. Если мы до- кажем, что даже в этом (наихудшем) случае характеристики системы останутся удовлетвори- тельными, то они будут не хуже и при любой другой случайной помехе. Корреляционная функция белого шума равна ) ( ) ( 0 τ δ τ S R X = . Действительно, преобразо- вание Фурье сразу дает ∫ ∞ ∞ − − = = 0 0 ) ( ) ( S d e S S j X τ τ δ ω ωτ Значения такого сигнала отстоящие по времени на любой, сколь угодно малый интервал, не- коррелированы. Это означает, что нет никакой зависимости между соседними, сколь угодно близко расположенными друг к другу, отсчетами такого случайного процесса. Корреляционная функция и спектральная плотность белого шума показаны на рисунках: Белый шум, как сигнал с бесконечной энергией, невозможно получить на практике. При моделировании его обычно заменяют на белый шум с ограниченной полосой, который имеет равномерный спектр в полосе частот от 0 ω − до 0 ω , и нулевой вне этой полосы: 5 Это название связано с белым светом, спектр которого содержит все частоты видимого спектра. 0 ω ) ( ω X S 0 ) ( S S X = ω 0 S 0 τ ) ( τ X R ) ( ) ( 0 τ δ τ S R X = 0 ω ) ( ω X S [ ] ) ( ) ( 2 ) ( 0 0 2 ω ω δ ω ω δ π ω − + + = A S X 0 ω 0 ω − 0 τ ) ( τ X R 2 / 2 A τ ω τ 0 2 cos 2 ) ( A R X = © К.Ю. Поляков, 2009 18 ⎪⎩ ⎪ ⎨ ⎧ > ≤ = 0 0 0 , 0 , ) ( ω ω ω ω ω S S X Средний квадрат такого сигнала равен π ω / 0 0 2 S x = , а не бесконечности. Корреляционную функцию можно найти с помощью обратного преобразования Фурье: τ ω τ ω π ω ω ω π τ 0 0 0 0 sin ) ( 2 1 ) ( ⋅ = = ∫ ∞ ∞ − S d S R X X Поскольку 0 sin = k π при любом целом k , корреляционная функция равна нулю при всех 0 ω π τ k = , где 0 ≠ k – любое целое число, не равное нулю. Это значит, что значения, взятые из такого сигнала в моменты времени K , 3 , 2 , , 0 0 0 0 ω π ω π ω π , (выборка с периодом 0 ω π ) будут некор- релированы. 0 ω ) ( ω X S 0 S 0 ω 0 ω − 0 τ ) ( τ X R τ ω τ ω π ω τ 0 0 0 0 sin ) ( ⋅ = S R X 0 ω π 0 2 ω π 0 3 ω π 0 ω π − © К.Ю. Поляков, 2009 19 3. Оценка и моделирование случайных процессов 3.1. Оценка корреляционной функции В прикладных задачах часто нужно определить корреляционную функцию и спектраль- ную плотность по экспериментальным данным. При этом мы можем наблюдать и анализиро- вать только «кусок» реализации на временном интервале от нуля до некоторого T , поэтому для невозможно использовать усреднение по ансамблю. Остается надеяться на то, что процесс эр- годический, и применять усреднение по времени. Пусть известна реализация случайного процесса ) (t x на интервале от 0 до T . Для оценки (приближенного вычисления) корреляционной функции при T << ≤ τ 0 (то есть при положи- тельных τ , достаточно малых по сравнению с T ) можно использовать формулу dt t x t x T R T X ) ( ) ( 1 ) ( ˆ 0 τ τ τ τ + − = ∫ − . (3) Обратите внимание, что время усреднения равно τ − T , а не T , потому что только интервал ] ; 0 [ τ − T содержит как t , так и τ + t . К сожалению, точно вычислить этот интеграл невозможно, потому что мы не знаем математическую формулу для ) (t x . В реальности обычно известны только значения этой функции (выборка) в моменты ∆ ∆ ∆ N , , 2 , , 0 K , где ∆ – интервал между измерениями. Тогда ) ( ˆ τ X R можно приближенно подсчитать только для ∆ ∆ ∆ = M , , 2 , , 0 K τ (где N M << ) по формуле ) ( ) ( 1 1 ) ( ˆ 0 ∆ + ∆ ⋅ ∆ + − = ∆ ∑ − = i k x k x i N i R i N k X , N M i << = , , 1 , 0 K , в которой интеграл заменен на сумму. С теоретической точки зрения математическое ожидание такой оценки (при усреднении по ансамблю) совпадает с истинной корреляционной функцией, то есть это – несмещенная оценка. 3.2. Оценка спектральной плотности 3.2.1. Использование оценки корреляционной функции Предположим, что мы исследуем эргодический процесс и знаем одну реализацию ) (t x на интервале от 0 до некоторого T . Выше было показано, что по этим данным можно построить оценку корреляционной функции. Если бы мы знали полностью непрерывную корреляционную функцию ) ( τ X R , для оценки спектральной плотности можно было бы использовать преобразо- вание Фурье (формулу (2)): ∫ ∞ = 0 cos ) ( ˆ 2 ) ( ˆ τ ωτ τ ω d R S X X В реальности известны лишь значения ) ( ˆ ∆ i R X в отдельных точках, поэтому последнюю фор- мулу нужно перевести в дискретный вид, заменив интеграл на конечную сумму: ∑ = ∆ ∆ ∆ = M i X X i i R S 0 cos ) ( ˆ 2 ) ( ˆ ω ω . (4) Этот метод оценки спектральной плотности называют методом Блэкмана-Тьюки. © К.Ю. Поляков, 2009 20 К сожалению, такой подход не всегда дает удовлетворительные результаты. Дело в том, что мы знаем только часть корреляционной функции, для значений τ от 0 до ∆ = M m τ . Эта не- полнота знаний может очень существенно влиять на результаты оценки спектра, вплоть до того, что вычисления по формуле (4) могут дать для некоторых частот отрицательные значения спек- тральной плотность. Этого не может быть в принципе, потому что мощность сигнала (и любой его составляющей) не может быть отрицательной. 3.2.2. Окна Чтобы исправить ситуацию, нужно как-то «сгладить» незнание корреляционной функции при больших τ и сделать оценку спектральной плотности более надежной. Для этого исполь- зуются так называемые «окна» – четные функции, на которые умножается корреляционная функция перед тем, как применить к ней преобразование Фурье. Одно из простейших «окон» – окно Хэмминга: ⎪⎩ ⎪ ⎨ ⎧ < ≥ + = m m m h w τ τ τ τ τ πτ τ , 0 , cos 46 , 0 54 , 0 ) ( На рисунке слева показано окно Хэмминга, а справа – исходная оценка корреляционной функ- ции ) ( ˆ τ X R и результат применения к ней окна Хэмминга ) ( ˆ ) ( τ τ X h R w (красная линия): Ясно видно, что применение этого окна (и других тоже) практически не изменяет форму корреляционной функции при малых τ , но сглаживает все выбросы при больших τ , которые, скорее всего, вызваны случайными ошибками. Для оценки спектральной плотности с учетом окна ) ( τ w применяют формулу, аналогич- ную (4): ∑ = ∆ ∆ ∆ ∆ = M i X X i i R i w S 0 cos ) ( ˆ ) ( 2 ) ( ˆ ω ω . (5) Не стоит печалиться по поводу того, что окно вносит дополнительное искажение. Так или иначе, «окно» используется всегда. Фактически, усекая корреляционную функцию, мы приме- няем прямоугольное окно: ⎩ ⎨ ⎧ < ≥ = m m r w τ τ τ τ τ , 0 , 1 ) ( На следующем рисунке показаны оценки спектра сигнала, полученные при использовании пря- моугольного окна ( ) ( ω X S , синяя линия) и окна Хэмминга ( ) ( ω h X S , красная линия). 0 τ ) ( τ X R m τ ) ( ˆ τ X R ) ( ˆ ) ( τ τ X h R w 0 τ ) ( τ h w 1 m τ m τ − © К.Ю. Поляков, 2009 21 Хорошо видно, что график ) ( ω X S заходит в отрицательную область, что невозможно с физиче- ской точки зрения. Применение окна Хэмминга позволило избавиться от этой проблемы и сгла- дить скачкообразные изменения оценки спектра. 3.2.3. Использование дискретного преобразования Фурье Главный недостаток классического метода оценки спектральной плотности (метода Блэк- мана-Тьюки) – большой объем вычислений. Гораздо меньше операций требуется при использо- вании прямого метода, основанного на использовании дискретного преобразования Фурье и со- временных вычислительных алгоритмах быстрого преобразования Фурье. При этом не нужно строить корреляционную функцию, а можно сразу найти спектральную плотность, обработав выборку значений исходного сигнала. В теории обработки аналоговых сигналов для перехода из временной области в частотную используется преобразование Фурье ∫ ∞ ∞ − − = dt e t f F t j ω ω ) ( ) ( Оно имеет смысл для любой детерминированной (неслучайной) функции ) (t f , которая абсо- лютно интегрируема, то есть интеграл от ее модуля на всей оси сходится: ∫ ∞ ∞ − ∞ < dt t f ) ( Для стационарного случайного процесса, не равного нулю, это условие никогда не будет вы- полняться, поэтому использовать преобразование Фурье в обычном смысле для анализа спектра случайных процессов нельзя. Однако если рассмотреть усеченный процесс ) (t x T , равный реализации ) (t x случайного процесса на интервале ] ; [ T T − и нулю вне этого интервала, для него можно найти преобразова- ние Фурье: ∫ ∫ − − ∞ ∞ − − = = T T t j t j T X dt e t x dt e t x F ω ω ω ) ( ) ( ) ( Квадрат модуля этой функции, деленный на ширину интервала T 2 , характеризует среднюю мощность сигнала на частоте ω . В пределе, при ∞ → T , мы должны получить спектральную плотность мощности. Так как ) (t x – это только одна реализация случайного процесса, в окон- чательной формуле нужно использовать усреднение по ансамблю (математическое ожидание): { } T F E S X T X 2 ) ( lim ) ( 2 ω ω ∞ → = 0 ω ) ( ω X S ) ( ω X S ) ( ω h X S © К.Ю. Поляков, 2009 22 При реальных измерениях мы знаем только одну реализацию случайного процесса ) (t x на интервале ] ; 0 [ T , поэтому усреднение по ансамблю чаще всего невозможно. Тогда для оценки спектральной плотности можно использовать предыдущую формулу без усреднения:: T F S X X 2 ) ( ) ( ˆ ω ω = , где ∫ − = T t j X dt e t x F 0 ) ( ) ( ω ω Теперь остается найти (приближенно) ) ( ω X F по дискретным измерениям процесса ) (t x Предположим, что известны его значения ) ( t k x x k ∆ ⋅ = при t k t ∆ ⋅ = для 1 ,..., 1 , 0 − = N k , так что интервал ] ; 0 [ T разделен на N подынтервалов шириной t ∆ (поэтому t N T ∆ ⋅ = ). Тогда интег- рирование можно приближено заменить суммой: ∑ ∑ − = ∆ ⋅ − − = ∆ ⋅ − ∆ = ∆ ⋅ ≈ 1 0 1 0 ) ( N k t k j k N k t k j k X e x t t e x F ω ω ω Для оценки спектра в теории обработки сигналов обычно используют сетку частот (в герцах) 1 , , 1 , 0 , − = ∆ ⋅ = N m t N m f m K с шагом t N f ∆ ⋅ = ∆ 1 . В теории управления принято строить спектры как функции угловой час- тоты (в радианах в секунду), которая получается из «обычной» частоты умножением на π 2 : 1 , , 1 , 0 , 2 − = ⋅ ∆ ⋅ = N m m t N m K π ω Для частоты m ω получаем m N k mk N j k m X X t e x t F ⋅ ∆ = ⋅ ∆ = ∑ − = − 1 0 2 ) ( π ω , (6) где через m X обозначена сумма, называемая дискретным преобразованием Фурье (ДПФ): ∑ − = − = 1 0 2 N k mk N j k m e x X π Заметим, что эта величина – комплексная, содержащая как вещественную, так и мнимую части. Легко подсчитать, что при расчете ДПФ по этим формулам для N частот количество опе- раций сложения и умножения будет пропорционально 2 N (обозначается ) ( 2 N O ). Это значит, что если N увеличивается, скажем, в 10 раз, то количество операций – примерно в 100 раз. Для больших N , особенно при анализе сигналов в реальном времени, такие расчеты выполняются недопустимо долго. Для быстрого вычисления ДПФ были разработаны специальные алгоритмы, которые на- зываются быстрым преобразованием Фурье (БПФ). Они позволили сократить количество опе- раций с ) ( 2 N O до ) log ( N N O . В функции |