Матлаб. К. Ю. Петрова введение в matlab учебное пособие
 Скачать 2.57 Mb.
|
|
minreal: >>mSYS=minreal(SYS) 4states removed. MATLAB сообщает, что порядок системы понижен на четыре и возвращает систему mSYS – минимальную реализацию пятого порядка. Заметим, что в данном случае из нуль-полюсного описания модели не удается усмотреть ни минимальной реализации, ни ее размерности. Более того, zpk-описания исходной и минимальной реализации совпадают: zpk(SYS), zpk(mSYS) Zero/pole/gain from input to output 1 Zero/pole/gain from input to output 2 Zero/pole/gain from input to output 3 1 #1: ----------------- (s+1) (s+2) (s+3) 1 #2: ----------------- (s+2) (s+3) (s+4) 1 #3: ----------------- (s+3) (s+4) (s+5) Команда minreal с двумя выходными параметрами [msys,U] = minreal(sys) дополнительно возвращает ортогональную матрицу преобразования U, что дает возможность найти каноническую декомпозицию Калмана исходной системы, используя формулы U*A*U',U*B,C*U'. Дискретные модели Описание дискретных моделей В MATLAB имеется возможность моделировать линейные системы с дискретным временем. Для описания таких систем используются разностные уравнения различных порядков, модели в пространстве состояний и дискретные передаточные функции. Линейное разностное уравнение n-го порядка имеет вид 0 1 1 1 1 0 1 1 1 1 k k k n n k k k n n k n u b u b u b y a y a y a y Здесь k – номер такта дискретного времени, u k – входной сигнал, у k – выходной сигнал, a i , b i – постоянные коэффициенты. Такими уравнениями, в частности, описываются цифровые фильтры. В ядре MATLAB имеется команда filter, позволяющая рассчитать выходной сигнал фильтра по известному входному сигналу, заданному массивом своих значений. Ее входными аргументами являются векторы а и b коэффициентов разностного уравнения, а также массив значений сигнала и: y=filter(b,a,u). С помощью четвертого входного аргумента можно задавать начальные условия фильтра. Гораздо больше возможностей по синтезу и моделированию цифровых фильтров предоставляет тулбокс SIGNAL, подробнее об этом говорится в конце данного раздела. Описание дискретной системы в пространстве состояний, как и в случае непрерывных систем, задается четверкой матриц: , 1 k k k k k k Du CX y Bu AX X 86 Здесь n k R X – вектор состояния; k k y u , – входной и выходной сигналы, A, B, C, D – постоянные матрицы. При исследовании линейных систем с дискретным временем и при решении разностных уравнений широко применяется математический аппарат z-преобразования, которое представляет собой дискретный аналог преобразования Лапласа. Дискретной передаточной функцией линейной системы с входом x(t) и выходом y(t) называется отношение z-преобразования выходного сигнала к z-преобразованию входного сигнала при нулевых начальных условиях Q(z) = Y(z)/U(z). Для систем с одним входом и одним выходом она представляет собой отношение двух полиномов ) ( 0 1 1 1 0 1 1 1 a z a z a z b z b z b z Q n n n n n Дискретная передаточная функция системы, заданной описанием в пространстве состояний, может быть найдена по формуле ) ( ) ( 1 D B A zE C z Q Дискретные модели во многом похожи на модели с непрерывным временем. Подробнее с ними можно ознакомиться по учебному пособию [4]. В табл. 4.3 приводится сравнение непрерывных и дискретных моделей. Таблица 4.3 Непрерывные Дискретные Вид уравнений Дифференциальные уравнения Разностные уравнения Описание в пространстве состояний Du CX y Bu AX X k k k k k k Du CX y Bu AX X , 1 Описание в виде передаточной функции , ) ( 0 1 1 1 0 1 1 1 a z a z a z a b z b z b z b z Q n n n n m m m m где z / 1 – оператор задержки , ) ( 0 1 1 1 0 1 1 1 a p a p a p a b p b p b p b p Q n n n n m m m m где p / 1 – оператор интегрирования Определение ПФ ) ( ) ( ) ( z U z Y z Q D B A pE C p U p Y p Q 1 ) ( ) ( ) ( ) ( Вид интегрального преобразования Преобразование Лапласа 0 ) ( ) ( dt e t f p F pt z -преобразование 0 ) ( ) ( n n z n f z F Свободное движение 0 ) ( X e t X At 0 X A X k k Компоненты свободного движения t i e , где i – простой корень характеристического уравнения; t m m j j j j e t k t k k ) ( 1 1 0 , где i – корень кратности m n i z , где i z – простой корень характеристического уравнения; n j m m j j j z n k n k k ) ( 1 1 0 , где j z – корень кратности m Тулбокс CONTROL позволяет моделировать линейные системы с непрерывным и дискретным временем. Дискретные модели в нем задаются при помощи тех же конструкторов, что 87 и непрерывные с той разницей, что последним параметром задается частота дискретизации ts (sampling time). Пример. Рассмотрим дискретную систему, описываемую уравнениями , 1 0 , 0 1 1 2 / 1 2 / 1 1 1 X y u X X k k k k при шаге дискретизации, равном 0,7. Создадим соответствующую ss-модель с дискретным временем: >> a=[-1 -1/2 ;1/2 -1]; b=[1;0]; c=[0 1]; d=0; ts=0.7; s1=ss(a,b,c,d, ts). На дисплей будут выведены матрицы: a = -1 -0.5 0.5 -1 b = 1 0 c = 0 1 d = 0 и сообщение: Sampling time: 0.7 Discrete-time model. Аналогично создаются tf и zpk-модели с дискретным временем, например: >> num=1; den=[1 1]; s2=tf(num,den, ts); >> z=[ ]; p=[-1 -2]; k=3; s3=zpk(z,p,k, ts) 1 1 2 z s 2 3 3 2 z z s Обращение ко всем полям дискретной модели осуществляется так же, как и у непрерывной. Время дискретизации хранится в поле ts : >> s3.ts ans = 0.7000 Моделирование дискретных систем Для моделирования дискретных систем используются те же команды (step, impulse, initial, lsim), что и в непрерывном случае. При получении реакции дискретной системы на входное воздействие с помощью команды lsim нужно помнить, что вектор временных отсчетов следует задавать, согласуясь с частотой дискретизации или же вовсе опускать. >> u=1./1:10; % входной сигнал вида 1 ½ 1/3 ….. 1/10. >> % Правильно: >> y=lsim(s1,u); y’ ans = 0 0 0.5000 0 0.8750 0.2500 0.9063 0.8750 0.6172 1.6719 >> % Неправильно: >> t=0:0.1:0.9; y1=lsim(s1,u,t) ??? Error using ==> rfinputs Time sample spacing must match sample time of discrete-time models. >> % Правильно: >> t=0:0.7:(0.7*9); y1=lsim(s1,u,t); Как было сказано, для получения реакции линейной дискретной системы на входной сигнал наряду с командой lsim тулбокса CONTROL можно использовать команду filter ядра MATLAB. Повторим предыдущий пример, используя эту команду: >> [n,d]=tfdata(s1); 88 >> yy=filter(n{1},d{1},u) yy = 0 0 0.5000 0 0.8750 0.2500 0.9063 0.8750 0.6172 1.6719 Видим, что результаты совпадают. При получении весовой функции дискретных систем вместо реакции на дельта-функцию рассматривают реакцию на единичный импульс вида . Для этого используется та же команда impulse. Реакцию на единичный скачок (переходную функцию) получают также при помощи команды step. Вместо функции freqs для дискретных систем используют функцию freqz. Например, для построения амплитудно-частотной характеристики дискретного звена с передаточной функцией ) 2 /( 1 ) ( z z Q нужно набрать: >> [h,w]=freqz(1,[1 2], 200); plot(w,abs(h)), grid on Любую непрерывную систему путем дискретизации времени можно приближенно преобразовать в дискретную. В Control Toolbox для преобразования непрерывной модели в дискретную и обратно существуют команды с2d (читается continuous to discrete) и d2c (discrete to continuous). Входным аргументом является исходная система. В команде с2d используется второй входной аргумент – время дискретизации. Допускается использование третьего входного аргумента, позволяющего выбрать один из шести возможных методов преобразования. Для этих же целей служат команды impinvar и bilinear (Signal Processing Toolbox). Ряд возможностей по работе с дискретными моделями и решению разностных уравнений содержит тулбокс SYMBOLIC. В нем имеются команды ztrans и iztrans для выполнения прямого и обратного z–преобразований. Например, в результате выполнения команд >>syms t, ztrans(t^2) получим ответ: ans =z*(z+1)/(z–1)^3, а в результате выполнения >>syms a t, ztrans(sin(a*t)), получим ответ: ans =z*sin(a)/(z^2-2*z*cos(a)+1). Тем самым мы нашли два табличных z-преобразования: 1 cos 2 sin sin , ) 1 ( ) 1 ( 2 3 2 a z z a z at z z z t Рассмотрим еще два примера дискретных моделей. Пример 1. Задача о банковском вкладе. Клиент открывает счет в банке, внося Врублей. Каждый год он добавляет по b рублей, годовой банковский процент равен k. Какая сумма окажется на его счете через n лет? Решение. Обозначая текущую величину вклада через y n , можно записать следующее разностное уравнение: y n+1 =y n (1+k)+b; y 0 =B. Применяя к нему z-преобразование ), B ) ( Y ( ); ( Y 1 z z y z y n n получаем алгебраическое уравнение: 1 ) ( Y ) 1 ( ) B ) ( Y ( z z b z k z z Выражаем из него Y(z): ) 1 )( 1 ( ) ( Y z k z b B Bz z z Выполним обратное z-преобразование с помощью команды iztrans: >> syms z b B k , Y=z*(B*z-B+b)/(z-k -1)/( z-1); y=iztrans(Y) y =(-b+(k+1)^n*B*k+(k+1)^n*b)/k 89 Таким образом, через n лет на банковском счете окажется сумма k b k k b y n n ) 1 )( B ( Возможности тулбоксов CONTROL и SYMBOLIC при работе с дискретными системами дополняют друг друга. Рассмотрим пример на их совместное применение. Пример 2. Дискретная система задана уравнениями 7 6 5 1 1 k k k k k k u y x x x y Требуется найти двумя способами (символьным и численным) передаточную функцию от входа u до выхода у и импульсную весовую функцию. Решение. В данном случае имеем следующие матрицы описания в пространстве состояний ]. 1 0 [ , 7 0 , 5 6 1 0 , C B A x y X Способ 1 (символьный). Дискретную передаточную функцию получаем по формуле , ) ( 1 B A zE C Q импульсную весовую функцию находим с помощью обратного z- преобразования. Проделаем эти выкладки в тулбоксе SYMBOLIC. Соответствующая программа имеет вид: >>syms z; >>А=[0 1;-6 5];В=[0; 7];С=[0 1]; %ввод матриц >>Q =С*inv(z*eye(2)-А))*В; %передаточная функция >>q=iztrans(Q); Q,q %весовая функция Q=7*z/(z^2-5*z+6),q=-7*2^n+7*3^n %результат Способ 2 (смешанный). Сначала используем команды тулбокса CONTROL: >> A=[0 1;-6 5];B=[0; 7];C=[0 1]; >>sd=ss(A,B,C,0,1); % Discrete-time model Sampling time: 1, >> s=tf(sd) Transfer function: 7 z ------------- z^2 - 5 z + 6 Переходим к символьному представлению: >> syms z; q=iztrans(7*z/(z^2-5*z+6)) q = -7*2^n+7*3^n В обоих случаях получаем одинаковый результат: , 6 5 7 ) ( 2 z z z z Q ). 2 3 ( 7 n n n q Расчет аналоговых и цифровых фильтров В пакете MATLAB предусмотрены средства для цифровой обработки сигналов и синтеза классических аналоговых и цифровых фильтров, таких как фильтры Баттерворта, Чебышёва и эллиптические фильтры. Аналоговые фильтры. В тулбоксе SIGNAL имеются специальные команды, позволяющие синтезировать аналоговые фильтры и осуществлять обработку сигналов с их помощью. В первую очередь здесь следует назвать команды buttap, cheb1ap, cheb2ap, ellipap, предназначенные для построения аналоговых фильтров-прототипов Баттерворта, двух разновидностей фильтров Чебышева и эллиптических фильтров. 90 Аналоговыми фильтрами-прототипаминазывают фильтры низких частот с единичной частотой среза. В дальнейшем они могут использоваться для построения других типов фильтров, например, высокочастотных, режекторных и т.д. Аналоговые низкочастотные прототипы фильтров реализуются с помощью следующих команд: buttap – прототип фильтра Баттерворта; cheb1ap – прототип фильтра Чебышева первого типа (равноволновые колебания в полосе пропускания); cheb2ap – прототип фильтра Чебышева второго типа (равноволновые колебания в полосе подавления); ellipap – прототип эллиптического фильтра. Передаточные функции фильтров Баттерворта и Чебышёва первого рода не имеют нулей: , ) ( 1 ) ( p A p Q n где ) ( p A n – некоторый полином n-го порядка. Квадраты их АЧХ ) ( ) ( * j Q j Q имеют вид n 2 1 1 и ) ( 1 1 2 2 n T соответственно, где T n – полином Чебышева порядка n. Полюсы фильтра Чебышева первого рода расположены в левой половине эллипса на s- плоскости. В пакете MATLAB аналоговый фильтр-прототип Чебышева первого рода рассчитывается с помощью команды [z, p, k] = cheb1ap(n, Rp), где n – порядок фильтра, Rp – уровень пульсаций в полосе пропускания (в децибелах). Пример. Синтезируем низкочастотный фильтр Чебышева восьмого порядка с единичной частотой среза и единичной величиной пульсаций в полосе пропускания: >> [z,p,k]=cheb1ap(8,1); Построим графики его амплитудно-частотной и фазочастотной характеристик: >>w=0:0.01:4; >>h=freqs(k*poly(z), poly(p),w); >>plot(w,abs(h)); >>sys=zpk(z,p,k); nyquist(sys) Оба графика показаны на рис. 5.3 0 0.5 1 1.5 2 2.5 3 0 0.2 0.4 0.6 0.8 1 w A -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 Nyquist Diagram Real Axis Im a g in a r y A x is Рис. 5.3 |