ОТУ_Поляков_Часть2_2009. Теория автоматического управления для чайников Часть ii. Управление при случайных
 Скачать 1.41 Mb.
|
|
fft среды M ATLAB используется модификация ал- горитма БПФ, предложенного Дж. Кули и Дж. Тьюки. Этот алгоритм наиболее эффективен, ес- ли число отсчетов N представляет собой степень двойки ( p N 2 = при целом 0 > p ). Заметим, что если это не так, всегда можно дополнить ряд нулями до ближайшей степени двойки. Казалось бы, формула (6) позволяет оценить спектр для всех частот вплоть до t N N N ∆ ⋅ − = − ) 1 ( 2 1 π ω . Однако нужно учесть, что для анализа мы используем только дискретные из- © К.Ю. Поляков, 2009 23 мерения с периодом t ∆ . Остальные значения непрерывного сигнала ) (t x (между моментами измерений) теряются, и с ними теряется информация о высокочастотных составляющих. Согласно теореме Котельникова-Шеннона, по дискретным измерениям с периодом t ∆ можно восстановить частотные свойства сигнала только до частоты t f ∆ = 2 1 max (или до соот- ветствующей угловой частоты t ∆ = π ω max , которая называется частотой Найквиста 6 ). Поэтому только оценка спектра на частотах 2 / 0 , , N ω ω K дает нам практически полезную информацию 7 Подведем итог. Для оценки спектра сигнала по N отсчетам ) 1 , , 1 , 0 ( − = N k x k K нужно выполнить следующие действия: 1) с помощью БПФ (функция fft в M ATLAB ) найти массив ) 1 , , 1 , 0 ( − = N m X m K ; 2) взяв первую половину этого массива, рассчитать соответствующие значения ) 2 / , , 1 , 0 ( ) ( N m X t F m m X K = ⋅ ∆ = ω для частот, не превышающих частоту Найквиста t N ∆ = = π ω ω 2 / max ; 3) для каждой частоты ) 2 / , , 1 , 0 ( N m m K = ω найти оценку спектральной плотности мощности по формуле: t N F T F S m X m X m X ∆ ⋅ = = 2 2 ) ( ) ( ) ( ˆ ω ω ω Для сглаживания спектральной плотности так же, как и в методе Блэкмана-Тьюки, ис- пользуются окна. Только теперь на весовую функцию умножается не оценка корреляционной функции, а сама реализация на интервале ] ; 0 [ T : Для этого случая окно Хэмминга на интервале ] ; 0 [ T принимает вид ⎪⎩ ⎪ ⎨ ⎧ > < ≤ ≤ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ − + = T t t T t T t T t w h , 0 , 0 0 , 2 2 cos 46 , 0 54 , 0 ) ( π Далее дискретное преобразования Фурье вычисляется для отсчетов взвешенной функции, то есть, вместо (6) получаем ∑ − = − ⋅ ∆ = 1 0 2 ) ( N k mk N j k k m X e w x t F π ω , где ) ( t k w w h k ∆ ⋅ = Использование окна для исходного сигнала приводит к уменьшению его энергии и, как следст- вие, к заниженным оценкам спектральной плотности. Чтобы скомпенсировать эти потери, весо- 6 Иначе говоря, синусоиду нужно измерять более 2 раз за период. 7 Можно доказать, что * k k N X X = − , где звездочка обозначает комплексно-сопряженную величину. 0 t ) ( ) ( t w t x h T 0 t ) (t x T 0 t ) ( τ h w T © К.Ю. Поляков, 2009 24 вая функция умножается на дополнительный коэффициент q , который часто выбирают из ус- ловия нормировки (сохранения энергии весовой функции окна, которая должна остаться такой же, как для прямоугольного окна): ∫ ∫ = = T T h T dt dt t w q 0 0 2 2 1 ) ( Несложно подсчитать, что для окна Хэмминга из этого условия следует 586 , 1 2 / 46 , 0 54 , 0 1 2 2 ≈ + = q 3.3. Прохождение случайных сигналов через линейные системы Существует два подхода к исследованию систем управления при случайных возмущениях: 1) вероятностный – на основе плотностей распределения вероятностей; 2) статистический – с помощью усредненных характеристик: математического ожидания, дисперсии, корреляционной функции, спектральной плотности. Применение вероятностного подхода, как правило, связано со значительными трудностя- ми. С одной стороны, они вызваны недостатком информации о плотностях распределения слу- чайных сигналов. С другой стороны, существующий математический аппарат достаточно сло- жен. Приведем только один важный факт: если входной сигнал имеет нормальное распределе- ние, то на выходе линейной системы будет также сигнал с нормальным распределением (ли- нейная система не «портит» нормальность). В прикладных задачах нас чаще всего интересует не плотность распределения вероятно- стей на выходе системы, а некоторые более осязаемые характеристики – среднее значение, дис- персия и т.д. Поэтому в подавляющем большинстве случаев используется статистический под- ход. Далее мы будем предполагать, что на вход линейной системы с известной передаточной функцией ) (s F действует стационарный случайный процесс с заданной спектральной плотно- стью ) ( ω X S Прежде всего, отметим, что при стационарном случайном входе выход ) (t y линейной стационарной системы (у которой все характеристики не зависят от времени) – тоже стацио- нарный случайный процесс. Для процесса ) (t y требуется найти • математическое ожидание y ; • дисперсию y D ; • корреляционную функцию ) ( τ Y R ; • спектральную плотность ) ( ω Y S Проще всего решается вопрос с математическим ожиданием: среднее значение выхода равно среднему значению входа, умноженному на статический коэффициент усиления системы (коэффициент усиления постоянного сигнала): ) ( lim , 0 s F k x k y s F F → = = Учитывая, что в линейных системах справедлив принцип суперпозиции (реакция на сумму двух сигналов равна сумме реакций на отдельные сигналы), далее мы для простоты будем рассмат- ривать только центрированные процессы, имеющие нулевые средние значения. © К.Ю. Поляков, 2009 25 Остальные характеристики удобнее определять с помощью спектральной плотности вы- хода. Вспомним, что спектральная плотность – это плотность распределения мощности сигнала по частотам. Сначала рассмотрим, как изменяется мощность гармонического сигнала t A t x 0 sin ) ( ω = , когда он проходит через линейную систему. Из классической теории автомати- ческого управления известно, что на выходе устанавливается гармонический процесс с той же частотой, но с другой амплитудой и фазой: ) sin( ) ( 0 ϕ ω + = t B t y , причем его амплитуда B определяется по частотной характеристике ) ( 0 ω j F : ) ( 0 ω j F A B ⋅ = Мощность гармонического сигнала (средний квадрат) пропорциональна квадрату амплитуды: ) ( ) ( ) ( 0 0 2 2 0 2 2 ω ω ω j F j F A j F A B − ⋅ ⋅ = ⋅ = В последнем равенстве использовано свойство комплексного числа ) ( 0 ω j F – квадрат его модуля равен произведению этого числа на комплексно-сопряженное, то есть на ) ( 0 ω j F − Таким образом, мы с помощью простых рассуждений вышли на очень важный результат: при гармоническом входе с частотой 0 ω мощность сигнала на выходе линейной системы равна мощности входного сигнала, умноженной на ) ( ) ( 0 0 ω ω j F j F − . Учитывая, что это свойство справедливо на всех частотах, и заменив слово «мощность» на «спектральную плотность», по- лучаем формулу, позволяющую сразу найти спектр процесса на выходе: ) ( ) ( ) ( ) ( ω ω ω ω j F j F S S X Y − ⋅ ⋅ = Соответствующая корреляционная функция ) ( τ Y R может быть найдена как обратное пре- образование Фурье от ) ( ω Y S . Для вычисление среднего квадрата процесса ) (t y нужно проин- тегрировать ) ( ω Y S : ∫ ∞ = 0 2 ) ( 1 ω ω π d S y Y Если процесс центрированный, средний квадрат совпадает с дисперсией, то есть 2 y D y = В общем случае для вычисления дисперсии нужно использовать равенство ( ) 2 2 y y D y − = Выделим один важный случай, когда входной сигнал – это единичный белый шум с посто- янной спектральной плотностью 1 ) ( = ω X S (белый шум единичной интенсивности). Тогда по- лучаем ) ( ) ( ) ( ω ω ω j F j F S Y − ⋅ = . (7) Таким образом, спектральная плотность выхода системы, на вход которой действует еди- ничный белый шум, равна квадрату ее амплитудной характеристики. Пусть передаточная функция линейной системы равна 1 1 ) ( + = s s F . Белый шум, проходя через такое звено, превращается в сигнал, имеющий спектральную плотность 1 1 ) ( ) ( ) ( 2 + = − = ω ω ω ω j F j F S Y , график которой показан на рисунке: © К.Ю. Поляков, 2009 26 Белый шум «содержит» все частоты, но они по-разному преобразуются. Постоянный сиг- нал (имеющий частоту 0 = ω ) передается на выход системы без изменений. На низких частотах искажения достаточно малые, а высокие частоты подавляются (фильтруются) системой. Это типичный фильтр низких частот (он пропускает низкочастотные сигналы и блокирует высо- кочастотные). Отметим очень важный факт: поскольку высокие частоты подавляются, отклоне- ния спектра входного сигнала от равномерного спектра белого шума в этой области не будут существенно влиять на спектр выхода. На этой идее основано компьютерное моделирование случайных процессов (см. далее). 3.4. Моделирование случайных сигналов К сожалению, анализ системы далеко не всегда можно выполнить теоретически. Это осо- бенно актуально для нелинейных систем. В этом случае единственным методом остается ими- тационное моделирование. Поэтому важно уметь моделировать случайные процессы, дейст- вующие на систему: возмущения (например, влияние ветра и волн на судно) и помехи измере- ния (погрешности измерительной системы). 3.4.1. Случайные числа Моделирование случайных процессов на цифровых компьютерах основано на использо- вании случайной последовательности чисел. Обеспечить подлинную случайность в программе практически очень сложно (иногда для этого используют шум звуковой карты, счетчик тактов процессора). Поэтому обычно применяют генераторы псевдослучайных чисел. В последова- тельности псевдослучайных чисел каждое следующее число 1 + k x рассчитывается по какой-то математической формуле на основе предыдущего k x (или нескольких предыдущих). Например, во многих системах программирования используется линейный конгруэнтный метод: m c ax x k k mod ) ( 1 + = + Здесь a , c и m – некоторые целые числа. По этой формуле вычисляются псевдослучайные числа, равномерно распределенные 8 на интервале от 0 до 1 − m . Для краткости будем далее на- зывать псевдослучайные числа просто случайными. Параметры a , c и m нужно выбрать так, чтобы полученная последовательность содер- жала как можно меньше закономерностей и как можно дольше не повторялась. Это отдельная математическая проблема, которая до сих пор не имеет однозначного решения. Например, в функции rand среды M ATLAB используется генератор с параметрами 9 2147483647 1 2 , 0 , 16807 7 31 5 = − = = = = m c a 8 Теоретически, конечно. Реальное распределение всегда немного неравномерное. 9 См. http://www.mathworks.com/company/newsletters/news_notes/pdf/Cleve.pdf 0 ω 1 1 ) ( 2 + = ω ω X S 1 ) ( ω X S © К.Ю. Поляков, 2009 27 Последовательность начинается с некоторого начального числа (англ. seed – семя, зерно). Чтобы получить в точности такую же последовательность (например, чтобы повторить полу- ченные ранее результаты моделирования), нужно начать последовательность с того же числа. Как показано в разд. 1.6, при сложении равномерно распределенных чисел распределение их суммы стремится к нормальному. Поэтому простейший способ получения нормального рас- пределения – сложить некоторое количество случайных чисел, полученных с помощью стан- дартного датчика. Часто используют суму 12 чисел, равномерно распределенных на интервале ] 1 ; 0 [ , потому что это сразу дает случайную величину с единичной дисперсией. В M ATLAB есть встроенная функция randn, которая генерирует последовательность слу- чайных чисел с нормальным распределением (с нулевым средним и единичной дисперсией) на основе более сложного алгоритма 10 3.4.2. Формирующий фильтр На практике обычно известна спектральная плотность ) ( ω X S и требуется смоделировать процесс, имеющий такую спектральную плотность. Вспомним, что спектральная плотность неотрицательна для любой частоты. Тогда функ- ция ) (s S X , полученная при подстановке 2 2 s − = ω в ) ( ω X S , неотрицательна на мнимой оси, то есть при ω j s = для всех ω . Можно доказать, что в этом случае ее можно представить в виде произведения ) ( ) ( ) ( s F s F s S X − = , то есть в форме (7). При этом всегда можно выбрать переда- точную функцию ) (s F так, чтобы она была устойчивой (не имела полюсов в правой полуплос- кости) и минимально-фазовой (не имела нулей в правой полуплоскости). Такой переход от ) (s S X к ) (s F называется факторизацией (англ. разложение на множители). Как следует из (7), если на вход звена с передаточной функцией ) (s F подать единичный белый шум, процесс на выходе будет иметь заданную спектральную плотность ) ( ω X S . Функ- ция ) (s F называется передаточной функцией формирующего фильтра. Проще всего моделировать процессы с дробно-рациональной спектральной плотностью. Например, одна из моделей морского волнения (подробности см. в главе 4) описывается спек- тральной плотностью 2 2 2 2 2 2 4 2 ) ( ) ( 2 4 ) ( β α ω β α ω αω ω + + − + = r D S , где r D – дисперсия волновой ординаты, α – коэффициент затухания и β – частота максимума спектра. Заменяя 2 ω на 2 s − , получаем 2 2 2 2 2 2 4 2 ) ( ) ( 2 4 ) ( β α β α α + + − − − = s s s D s S r Очевидно, что формирующий фильтр будет иметь передаточную функцию вида 0 1 2 ) ( δ δ γ + + = s s s s F , так что 10 См. http://www.mathworks.com/company/newsletters/news_notes/clevescorner/spring01_cleve.html © К.Ю. Поляков, 2009 28 2 0 2 0 2 1 4 2 2 ) 2 ( ) ( ) ( δ δ δ γ + − − − = − s s s s F s F Приравнивая коэффициенты числителя и знаменателя ) (s S и ) ( ) ( s F s F − , сразу находим: α γ r D 2 = , 2 2 0 β α δ + = , α δ β α δ 2 ) ( 2 0 2 2 1 = + − = В более сложных случаях факторизация выполняется с помощью численных методов. Нужно разложить на простейшие сомножители числитель и знаменатель ) (s S и включить в ) (s F только те множители, корни которых находятся в левой полуплоскости (в числителе до- пускаются корни на мнимой оси). Посмотрим, как работает этот метод на том же примере. Чис- литель спектральной плотности ) (s S имеет двойной корень в точке 0 = s , его можно предста- вить в симметричном виде ) ( 2 2 4 2 s D s D s D r r r − ⋅ ⋅ ⋅ = − α α α Знаменатель имеет корни в точках β α j ± и β α j ± − . Учитывая, что 0 > α , определяем, что корни β α j ± − расположены в левой полуплоскости. Поэтому знаменатель ) (s F представляет собой произведение 2 2 2 2 ) )( ( β α α β α β α + + + = − + + + s s j s j s . Таким образом, ) ( ) ( ) ( ) ( 2 4 ) ( 2 2 2 2 2 2 4 2 s F s F s s s D s S r − = + + − − − = β α β α α , где 2 2 2 2 2 ) ( β α α α + + + ⋅ = s s s D s F r Такой ж результат был получен ранее методом неопределенных коэффициентов. Итак, формирующий фильтр мы построили. Теперь остается один очень практический во- прос: как получить белый шум, который, как известно, является сигналом с бесконечной энер- гией? Вспомним, что белый шум – это только вспомогательный сигнал, который, проходя через систему с передаточной функцией ) (s F , генерирует сигнал с заданной спектральной плотно- стью. Оказывается, можно заменить его на другой сигнал (который просто получить на компь- ютере), и при этом спектральная плотность выхода оказывается достаточно близка к заданной. 3.4.3. Случайный ступенчатый сигнал вместо белого шума Как мы видели, на компьютере несложно получить последовательность случайных (точ- нее – псевдослучайных) чисел с равномерным или нормальным распределением. По этим чис- лам можно построить ступенчатый сигнал, фиксируя каждое значение в течение некоторого времени k τ : Теоретически эти числа некоррелированы 11 ; при этом можно показать, что корреляцион- ная функция ) ( τ X R ступенчатого сигнала – треугольная (см. рисунок справа). При 0 = τ она равна дисперсии D последовательности случайных чисел, а при k τ τ > обращается в нуль (по- 11 Строго говоря, это не совсем так. Из-за неидеальности компьютерного датчика псевдослучайных чисел эти зна- чения будут коррелированны и корреляционная функция будет немного отличаться от треугольной. 0 t ) (t x k τ 0 τ ) ( τ X R D k k X D R τ τ τ τ τ < ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ − = , 1 ) ( k τ k τ − © К.Ю. Поляков, 2009 29 тому что моменты времени t и τ + t находятся на разных интервалах и, следовательно, соответ- ствующие значения некоррелированы). Число k τ называют интервалом корреляции – так назы- вается интервал, после которого можно считать корреляционную функцию (примерно) равной нулю. Взяв преобразование Фурье от корреляционной функции ⎪ ⎩ ⎪ ⎨ ⎧ ≥ < ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ − = k k k X D R τ τ τ τ τ τ τ , 0 , 1 ) ( получаем спектральную плотность k k X D S τ ω ωτ ω 2 ) cos 1 ( 2 ) ( − = Вычисляя предел этой функции при 0 → ω , находим k X D S τ ω ω = → ) ( lim 0 , так что при выборе k D τ / 1 = это значение равно 1 (как у белого шума). Заметим, что 0 ) ( = ω X S , когда 1 cos = k ωτ , то есть m k π ωτ 2 = при любом целом m . Форма спектральной плотности показана на рисунке ниже (здесь и далее принимается k D τ / 1 = ): Конечно, это далеко не белый шум, у которого спектральная плотность должна быть по- стоянной на всех частотах. Тем не менее, при уменьшении интервала корреляции k τ «колокол» расширяется, и для низких частот можно считать, что 1 ) ( ≈ ω X S . В пределе при 0 → k τ спектр стремится к равномерному спектру единичного белого шума. Далее будет показано, что при грамотном выборе k τ такой сигнал можно использовать в качестве источника вместо белого шума. Для примера предположим, что нужно получить сигнал со спектральной плотностью 1 1 ) ( 2 + = ω ω Y S , то есть формирующий фильтр имеет передаточную функцию 1 1 ) ( + = s s F . В ка- честве входного сигнала для этого звена будем использовать описанный выше ступенчатый сигнал при 5 , 0 = k τ с. На рисунке приведены графики спектральной плотности ступенчатого сигнала (синяя линия), желаемой спектральной плотности (сплошная зеленая линия) и фактиче- ской спектральной плотности выхода (штриховая линия). 0 ω 2 2 ) cos 1 ( 2 ) ( k k X S τ ω ωτ ω − = k τ π 2 − 1 k τ π 2 k k τ τ < 2 ) ( ω X S © К.Ю. Поляков, 2009 30 По графику видно, что в существенной полосе частот (где частотная характеристика звена ненулевая) спектр входного сигнала существенно неравномерный, поэтому желаемый и факти- ческий спектры на выходе системы немного различаются в области высоких частот. Приняв c k 01 , 0 = τ , имеем совершенно другую картину: Спектр входного сигнала в интересующей нас области практически равномерный, в спектр реального выхода практически точно совпадает с заданным. Очевидно, что при выборе k τ нужно учитывать частотные свойства формирующего фильтра, точнее, полосу частот, где его частотная характеристика достаточно отличается от ну- ля. Для этого используют понятие полосы пропускания системы b ω – так называется частота, для которой амплитудная частотная характеристика уменьшается на 3 дБ (децибела) в сравне- нии с максимальным значением (составляет примерно 0,708 от максимума). Разработчики M ATLAB рекомендуют при моделировании использовать значение b k ω π τ 2 100 1 ⋅ = В нашем случае амплитудной частотная характеристика (апериодического звена) имеет вид 1 1 ) ( 2 + = ω ω j F , ее максимум равен 1 (при 0 = ω ), поэтому полоса пропускания опреде- ляется равенством 708 , 0 1 1 2 = + b ω . Отсюда следует 998 , 0 1 708 , 0 1 2 ≈ − = b ω , так что 063 , 0 ≈ k τ 1 ) ( ω S c k 01 , 0 = τ ω 1 ) ( ω S c k 5 , 0 = τ ω |