Матлаб. К. Ю. Петрова введение в matlab учебное пособие
 Скачать 2.57 Mb.
|
|
dsolve тулбокса SYMBOLIC. Чтобы получить функцию h(t) для нашего примера, нужно набрать: >>syms a b; h=dsolve (‘Dx+a*x=b, x(0)=0’). На дисплее появится ответ h=b/a-exp(-a*t)*b/a, совпадающий с полученным выше. В библиотеке CONTROL весовая и переходная характеристики получаются путем прямого численного моделирования. Соответствующие команды называются impulse и step,это сокращения от impulse function (импульсная функция) и step function (переходная функция). Для выполнения этих операций во всех случаях необходимо предварительно ввести исходную информацию о системе в виде tf-модели или ss-модели, а также сформировать массив равноотстоящих моментов времени t, задающий временной интервал моделирования. Существует несколько модификаций команды impulse. Простейшая из них имеет вид impulse(sys), ее результатом является график весовой функции. Если набрать impulse(sys,10), то график будет построен на интервале 10 0 t (в предыдущем случае MATLAB сам определял длительность интервала). Более содержательные варианты этой команды получаются, если использовать выходные параметры (их число можно задавать от одного до трех). Наиболее полный вариант имеет вид [y,t,X]=impulse(sys,t), он предполагает, что sys – это ss-модель. Здесь t – массив точек времени, который нужно сформировать заранее (например, t=0:.1:10 ), y – выходной сигнал, Х – вектор состояния. Если моделируется система второго порядка на указанном интервале времени, то массив Х будет содержать два столбца по 101 числу в каждом. Первый столбец – это отсчеты функции х 1 (t), второй – отсчеты функции х 2 (t). Столбец у будет содержать 101 значение выходного сигнала у(t). Для построения графиков этих сигналов нужно использовать команды plot(t,y), plot(t,X), plot(t,X(:,1)) . В первом случае будет выведен график функции у(t), во втором – графики обеих переменных х 1 (t), х 2 (t), в третьем – график одной переменной х 1 (t). Команда plot(X(:,1), X(:,2)) построит траекторию на фазовой плоскости ). ( 1 2 x f x Функция step обеспечивает получение переходной функции модели, т.е. реакции на входной сигнал в виде единичной ступеньки. Она имеет те же модификации, что и impulse: >>step(sys), step(sys,T), y=step(sys,t), [y,t,X]=step(sys,t). Здесь, как и раньше, в качестве второго входного аргумента можно указывать либо число T (последний момент времени), либо массив t (все точки временного интервала). Аналогичный синтаксис имеет и команда initial (от initial condition – начальные условия). Она позволяет моделировать свободное движение системы, заданной своим матричным описанием, из начальных условий Х 0 (входной сигнал при этом не подается и=0). Просто initial(sys,X 0 ) строит график выходного сигнала, в более полном варианте используются три выходных параметра [y,t,X]=initial(sys,X0,t). 30 Основная команда, применяемая для моделирования линейных систем – это команда lsim (от linear simulation – линейное моделирование). Она обеспечивает получение реакции модели на произвольный входной сигнал и(t), представленный массивом своих отсчетов. Простейшая модификация этой команды lsim(sys,u,t) выводит график выходного сигнала системы. Естественно, предварительно надо сформировать массив времени t, входной сигнал и ввести описание системы sys. Например, задав: >>t=0:1:10; u=sin(t); sys=tf(1,[1 1]; lsim(sys,u,t), получим реакцию апериодического звена с передаточной функцией 1 1 p на входной сигнал u=sint при 10 0 t . Другие модификации этой команды могут использовать выходные параметры и начальные условия: >>y=lsim(sys,u,t), [y,t,X]=lsim(sys,u,t), [y,t,X}=lsim(sys,u,t,X0). Последний вариант требует задания системы sys в виде ss-модели. Построение графиков производится с помощью команды plot, например plot(t,u,t,y). В пакете MATLAB имеется еще ряд возможностей для расчета отклика линейных систем на известные входные воздействия. К ним относятся: – использование матричной экспоненты для получения весовой и переходной функций (команда expm); – вычисление интеграла свертки входного сигнала и весовой функции системы (команда trapz); – аналитическое решение дифференциальных уравнений (команда dsolve тулбокса SYMBOLIC); – использование прямого и обратного преобразований Лапласа (команды laplace и ilaplace тулбокса SYMBOLIC); – структурное моделирование в SIMULINK. Часть из них будет рассмотрена в следующих разделах. 2.3 . Частотные характеристики В библиотеке CONTROL имеются функции для исследования систем в частотной области. К основным частотным характеристикам относятся амплитудно-частотная характеристика (АЧХ), фазо-частотная характеристика (ФЧХ) и амплитудно-фазовая характеристика (АФХ), называемая также диаграммой Найквиста. Все они могут быть получены из передаточной функции системы Q(p) после формальной подстановки i p , где , 1 i – вещественная переменная, изменяющаяся в пределах от нуля до бесконечности. С физической точки зрения – это частота синусоидального сигнала, подаваемого на вход системы, а Q(i ) – Фурье изображение соответствующего выходного сигнала. Число ) ( ) ( ) ( ib a i Q – комплексное, оно может быть изображено на комплексной плоскости с помощью вектора с координатами (a, b) (рис. 2.3). Длина этого вектора определяется формулой , 2 2 b a A а угол φ – соотношением a b tg Если изменять частоту от нуля до бесконечности, конец вектора опишет некоторую траекторию на комплексной плоскости. Она называется Im Re b a Q(i ) A φ Рис. 2.3 31 амплитудно-фазовой характеристикой системы или годографом Найквиста (в последнем случае берут ). Зависимость длины вектора от частоты А( ) называется амплитудно-частотной характеристикой, а зависимость φ( ) – фазо-частотной характеристикой. Графики АЧХ и АФХ часто рисуют в логарифмическом или полулогарифмическом масштабе. Логарифмическая амплитудная характеристика (ЛАХ) описывается формулой ) ( log 20 10 A и изображается в логарифмическом масштабе. ЛАХ и АФХ, построенные в логарифмическом масштабе, называются в зарубежной литературе диаграммами Боде. Опишем процедуру получения этих характеристик в MATLAB. Будем считать, что система задана своей tf-моделью и сформирован вектор частот W, например: >>sys=tf([1 1], [1 1 1]); W=0: .1: 10; Для получения частотного отклика (frequency response) Q(i ) используется команда freqresp. В нашем примере передаточная функция и частотный отклик имеют вид 1 1 ) ( , 1 1 ) ( 2 2 i i i Q p p p p Q По команде H=freqresp(sys,W) будет вычислен массив комплексных чисел, содержащих 101 значение функции Q(i ). Команда plot(H(:)) построит зависимость мнимой части этих чисел от вещественной, т.е. график АФХ (рис. 2.4). 0 0.2 0.4 0.6 0.8 1 1.2 1.4 -1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 Рис. 2.4 32 0 2 4 6 8 10 0 0.5 1 1.5 0 2 4 6 8 10 -1.6 -1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 Рис. 2.5 Рис. 2.6 Для построения графиков АЧХ и ФЧХ надо найти модуль и аргумент этих чисел, для чего служат команды abs и angle: >>A=abs(H(:)); f=angle(H(:)); plot(W,A), plot((W,f). Результат показан на рис. 2.5, 2.6, из них видно, что с ростом частоты АЧХ стремится к нулю, а ФЧХ – к величине -π /2. Более короткий путь построения частотных характеристик – использование команд nyquist и bode. Команда nyquist(sys) построит диаграмму Найквиста, т.е. АФХ (при этом рассматриваются как положительные, так и отрицательные частоты). По команде bode(sys) будут построены графики АЧХ и ФЧХ в логарифмическом масштабе. Еще один способ построения графиков частотных характеристик, а также импульсной и переходной функций – это использование специального графопостроителя тулбокса CONTROL, который вызывается командой ltiview. 2.4. Анализ линейных систем К классическим задачам исследования линейных систем относятся отыскание нулей и полюсов передаточной функции, а также анализ устойчивости, ограниченности, управляемости и наблюдаемости. Коротко охарактеризуем команды библиотеки CONTROL, предназначенные для решения этих задач (табл. 1). Таблица 1 pole eig ctrb gram dcgain zero (tzero) pzmap obsv minreal roots Наиболее известный способ проверки устойчивости связан с вычислением полюсов системы, т.е. корней ее характеристического уравнения. Критерий устойчивости. Для того чтобы линейная система была устойчивой, необходимо и достаточно, чтобы все корни характеристического уравнения лежали в левой полуплоскости. Пример (анализ устойчивости и ограниченности). Дана неоднородная система дифференциальных уравнений 2 , 4 2 b y x y a x y x 33 Требуется проанализировать ее устойчивость и выяснить, при каких постоянных a и b все решения системы ограничены. Решение. Запишем систему в матричной форме B AX X , где b a B A , 1 2 2 4 и найдем корни ее характеристического уравнения, т.е. собственные числа матрицы А: >> A=[-4 2;2 -1]; eig(A) ans = -5 0 Одно из собственных чисел отрицательно, а другое лежит на мнимой оси, следовательно, однородная система находится на границе устойчивости. Для ответа на вопрос об ограниченности неоднородной системы найдем ее аналитическое решение с помощью команды dsolve : >> syms a b;s=dsolve('Dx=2*y-4*x+a,Dy=2*x-y+b') s = x: [1x1 sym] y: [1x1 sym] >> x=s.x x =2/5*exp(-5*t)*C1-1/10*b+1/5*a+2/5*t*b+1/5*t*a+1/2*C2 >> y= s.y y =-1/5*exp(-5*t)*C1+4/5*t*b+2/5*t*a+C2 Приведем подобные члены и перейдем к обычной нотации: 2 5 1 2 5 1 2 ) 2 ( 2 2 ) 2 ( c e c t b a y c a e c t b a x t t Решение будет ограниченным, если a=-2b, тогда 2 , 2 2 5 1 2 5 1 c e c y c a e c x t t Нули и полюсы системы, заданной передаточной функцией – это просто корни i z и i p полиномов, стоящих в числителе и знаменателе. Поэтому для вычисления вектора нулей z и вектора полюсов p передаточной функции Q(p) = num/den, могут использоваться команды z = roots(num); p = roots(den). Если система sys задана как tf-модель или ss-модель, то используются команды p=pole(sys), z=zero(sys). Для нахождения полюсов допустимо также использование команды p=eig(sys), что эквивалентно команде p=eig(sys.a), т.е. вычислению собственных чисел матрицы А. Функция tzero (от transfer zeros – передаточные нули) позволяет находить нули системы по матрицам описания в пространстве состояний z=tzero(A,B,C,D). Функция pzmap предназначена для одновременного вычисления нулей и полюсов. Если набрать [p, z]=pzmap(sys), то будут выведены столбцы p и z полюсов и нулей, а просто pzmap(sys) показывает расположение нулей и полюсов на комплексной плоскости (на графике нули изображаются ноликами, а полюсы – крестиками). При анализе управляемости и наблюдаемости линейных систем используются матрицы управляемости R и наблюдаемости D, построенные на основе матриц A, B , C описания в пространстве состояний: CA CA C D ], B A ..., AB, , [B R 1 1 n n Приведем соответствующие критерии. 34 Критерий управляемости. Для того чтобы система была управляемой, необходимо и достаточно, чтобы матрица управляемости имела полный ранг: rankR = n. Критерий наблюдаемости. Для того чтобы система была наблюдаемой необходимо и достаточно, чтобы матрица наблюдаемости имела полный ранг: rankD = n. Критерий минимальности. Для того чтобы система была минимальной необходимо и достаточно, чтобы обе матрицы R и D имели полный ранг: rankR = n, rankD = n. Формирование матриц управляемости и наблюдаемости производится с помощью команд ctrb и obsv (от controlability и observability). Их синтаксис одинаков: R=ctrb(sys), D=obsv(sys). В качестве аргументов можно использовать непосредственно матрицы А, В, С описания в пространстве состояний, например R=ctrb(A, B), D=obsv(A, C) или R=ctrb(sys.a, sys.b), D=obsv(sys.a, sys.c). Для вычисления ранга этих матриц используется команда rank. Другой способ проверки управляемости и наблюдаемости опирается на вычисление рангов так называемых грамианов управляемости и наблюдаемости. Для их нахождения служит команда gram, дополнительные сведения о ней можно найти в разд. 4.4, 4.5. Если система неуправляема или ненаблюдаема, то ее порядок может быть понижен путем удаления неуправляемых и ненаблюдаемых подсистем. У SISO-систем это эквивалентно сокращению совпадающих нулей и полюсов передаточной функции. Для этой цели используется команда minreal (от minimal realization). Ее можно использовать с одним и двумя входными аргументами: sys1=minreal(sys) и sys1=minreal(sys, eps). Второй аргумент позволяет указывать допуск ε, задающий степень близости сокращаемых нулей и полюсов. Пример (анализ минимальности). Рассмотрим систему с передаточной функцией 2 3 1 ) ( 2 p p p p Q Она имеет один нуль z 1 = – 1 и два полюса р 1 = – 1, р 2 = – 2. Система неминимальна, поскольку р 1 = z 1 . Сформируем tf-модель этой системы и найдем ее нули и полюсы: >>s=tf([1 1],[1 3 2]) >>zero(s) >>pole(s) >>eig(s) >> [p, z]=pzmap(s) s + 1 ------------- s 2 + 3 s + 2 -1 -2 -1 -2 -1 p = – 2 z = – 1 – 1 Расположение нулей и полюсов на комплексной плоскости можно получить, набрав команду pzmap(s) без выходного аргумента (на рис. 2.7, полюсы помечены крестиками, нуль – ноликом). 35 Оба полюса отрицательны, следовательно, система устойчива. Наличие диполя (совпадающего нуля и полюса) в точке (-1, 0) говорит о неминимальности системы. Перейдем к ss-модели и проанализируем ее управляемость и наблюдаемость, используя команды ctrb, obsv и rank: Рис. 2.7 >>s1 = ss(s) >>R=ctrb(s1) >>D=obsv(s1) >>rank(R) >>rankD a = -3 -2 1 0 b = 2 0 c = 0.5 0.5 R = 2 -6 0 2 D = 0,5 0,5 -1 -1 2 1 Анализируя ранги матриц R и D, заключаем, что система s1 управляема, но ненаблюдаема, следовательно, ее порядок может быть понижен. Найдем передаточную функцию минимальной реализации: q=minreal(s), получаем ответ q 2 1 s К тому же результату приходим, сокращая числитель и знаменатель исходной передаточной функции ) ( p Q на общий множитель р+1. Заметим, что статический коэффициент усиления при переходе к минимальной реализации не изменился: >>K=dcgain(s) >>K=dcgain(q) K=0.5 k=0.5 Весовая и переходная функции также остаются прежними. В этом можно убедиться с помощью команд impulse(s, q), step(s, q), по которым будут построены графики указанных функций для обеих систем. Команды ctrb и obsv можно использовать при работе с символьными выражениями, например, когда часть элементов матриц A, b, c заданы в буквенном виде. Пример (анализ управляемости и наблюдаемости системы третьего порядка). Объект управления задан описанием в пространстве состояний X = AX + bu, y = cX, 1 1 1 , 1 , 1 0 0 0 0 1 A 3 2 1 2 1 c a a b a a a , где X R 3 – вектор состояний, u, y – входной и выходной сигналы. Требуется проанализировать его управляемость и наблюдаемость. Решение. Вводим исходные символьные матрицы >> syms a1 a2 a3 real >> A=[-a1 1 0;0 -a2 0;0 1 -a1], b=[a2;a3;1]; с=[1 1 1]; -2 -1.5 -1 -0.5 0 -1 -0.5 0 0.5 1 Pole-Zero Map Real Axis Im a g in a r y A x is 36 Формируем матрицы управляемости и наблюдаемости R=[b, Ab, A 2 b], D=[c T , (cA) T , (cA 2 ) T ] T : >> R=ctrb(A, b) >> D=obsv(A, c), [ a2, -a1*a2+a3, -a1*(-a1*a2+a3)-a2*a3] [a3, -a2*a3, a2^2*a3 ] [1, a3-a1, -a2*a3-a1*(a3-a1)] [ 1, 1, 1 ] [ -a1, 2-a2, -a1] [a1^2, -2*a1-(2-a2)*a2, a1^2] Вычисляем определители: det(R) =0, det(D) =0. Оба определителя равны нулю, следовательно, система неуправляема и ненаблюдаема. Вырожденность матрицы наблюдаемости очевидна (ее первый и третий столбцы совпадают). Вырожденность матрицы управляемости «на глаз» обнаружить значительно сложнее. |