Лабораторная работа, система сглаживания. Лабораторная работа 5. Исследование разомкнутых ли нейных систем Цель работы
 Скачать 477.37 Kb.
|
|
Лабораторная работа №5. Исследование разомкнутых ли- нейных систем Цель работы Освоение методов анализа линейной непрерывной системы с помощью среды Scilab. Задачи работы 1. Научиться вводить описание линейной системы в среде Scilab.. 2. Получить навыки построения эквивалентных моделей в пространстве со- стояний и в форме передаточной функции. 3. Научиться строить средствами Scilab импульсную и переходную характе- ристики, карту расположения нулей и полюсов, амплитудно-частотную, фазо-частотную и амплитудно-фазовую характеристики. 4. Научиться строить график отклика линейной системы на произвольное входное воздействие. Краткие теоретические сведения Способы определения линейной системы Аналитически линейная система может быть задана либо в пространстве состояний (в виде уравнений состояния и выхода), либо передаточной функци- ей. Кроме аналитических существуют и графические формы (например, карта нулей и полюсов — полюсно-нулевое представление). Уравнения состояния в нормальной форме, как известно, имеют вид , (1) ⎧ = + ⎨ = + ⎩ X AX BU Y CX DU где X— вектор переменных состояния, U— вектор внешних воздействий (входных сигналов), Y— вектор выходных величин (сигналов). Передаточная функция W(p) (или W(s)) системы с одним входом и одним выходом представляет собой рациональную дробь, числитель и знаменатель которой представлены полиномами дифференциального оператора dt d p / = или комплексной частоты ω + α = j s . (Хотя смысл переменных p и s различный, сами полиномы выглядят абсолютно одинаково.) K( ) W( ) D( ) p p p = Полином K(p) (или K(s)) называют динамическим коэффициентом пере- дачи, а полином D(p) (D(s)) — характеристическим полиномом. Например, если линейная система описывается дифференциальным урав- нением вида , (2) 3 2 2 y y y x + + = + x где x(t) — входной сигнал, y(t) — выходной, то используя дифференциальный оператор dt d p / = , уравнение можно переписать в виде 2 3 2 2 p y py y px x + + = + или ( ) ( ) 2 3 2 1 2 1 p p y p + + = + x . Отсюда 2 ( ) 2 1 ( ) 3 2 1 y t p x t p p + = + + (3) Применив к тому же уравнению (2) преобразование Лапласа [ ] 0 ( ) ( ) ( ) st L f t F s f t e dt ∞ − = = ∫ в предположении нулевых начальных условий, получим ( ) ( ) 2 3 2 1 ( ) 2 1 ( ) s s Y s s X s + + = + . Отсюда 2 ( ) 2 1 W( ) ( ) 3 2 1 Y s s s X s s s + = = + + (4) здесь X(s) и Y(s) — изображения по Лапласу входного и выходного сигналов. Связь между передаточной функцией и математической моделью в про- странстве состояний задаётся известной формулой (5) ( ) 1 W( ) s s − = − + C 1 A B D Если входных и (или) выходных сигналов несколько, то передаточная функция будет представлять собой матрицу соответствующего размера, со- стоящую из рациональных дробей вида K ( ) W ( ) D( ) ij ij s s s = , где i — номер выхода, а j — номер входа. Характеристический полином ( ) D( ) det s s = − 1 A от количества входов и выходов не зависит и определяется только внутренними свойствами системы. Определить линейную систему в Scilab можно функцией syslin(), ко- торая возвращает указатель на структуру определённого формата и может вы- зываться несколькими способами. [sl]=syslin(dom,A,B,C [,D [,x0] ]) [sl]=syslin(dom,Num,Den) [sl]=syslin(dom,H) Параметр dom (domain) задаёт временной аспект управления и может принимать следующие значения: dom='c' для систем непрерывного (continuous) управления, 2 dom='d' для дискретных (discrete) систем управления, причём задание веще- ственного числа dom=n (как вариант) определяет дискретную систему с перио- дом дискретизации в n секунд. dom=[] (не задан), если временной аспект не определён. A,B,C и D определяют матричные коэффициенты математической моде- ли в пространстве состояний, а x0 — начальные условия (значения переменных состояния в начальный момент времени). Num, Den — полиномы (или матрицы полиномов) числителя и знамена- теля передаточной функции. H — рациональная дробь (или матрица рациональных дробей) передаточ- ной функции. Выполните следующий пример. A=[0,1;1,-1] B=[1;1] C=[1,1] D=1 S1=syslin('c',A,B,C,D) Изучите структуру, на которую указывает S1: S1(1) S1(2) //или S1("A") S1(3) //или S1("B") S1(4) //или S1("C") S1(5) //или S1("D") S1(6) //или S1("X0") S1(7) //или S1("dt") В графическом окне 1 постройте диаграмму Боде (АЧХ и ФЧХ) линейной системы S1 в диапазоне частот (0.01…10)Гц: scf(1); clf(1); bode(S1, 0.01, 10) Теперь определите ту же систему через передаточную функцию: s=poly(0,'s'); E=eye(); H=C*((E*s-A)^-1)*B+D S2=syslin('c',H) Изучите структуру, на которую указывает S2: S2(1) S2(2) //или S2("num") S2(3) //или S2("den") S2(4) //или S2("dt") Постройте диаграмму Боде системы S2 в графическом окне 2: scf(2); 3 clf(2); bode(S2, 0.01, 10) Определите линейную систему S3 через полиномы числителя и знамена- теля передаточной функции и постройте диаграмму Боде в графическом окне 3: K3=numer(H) D3=denom(H) S3=syslin('c',K3,D3) scf(3); clf(3); bode(S3, 0.01, 10) Сравните графические окна 1, 2 и 3. Если Вы были внимательны, то заметили, что ссылки S1 и S2 указывают на разные структуры, хотя линейные системы, которые они определяют, абсо- лютно идентичны. Структура S1 описывает систему в пространстве состояний, а S2 — через передаточную функцию. В Scilab имеется возможность преобра- зования друг в друга этих структур. Например, чтобы перейти от системы, за- данной в пространстве состояний к системе, заданной передаточной функцией, можно воспользоваться встроенной функцией ss2tf() (state-space to transfer function). Для обратного преобразования используется функция tf2ss(). Правда, в отличие от перехода ss2tf, который однозначен, обратный переход tf2ss неоднозначен, т.е. одной и той же передаточной функции могут соот- ветствовать множество уравнений состояния и выхода. В этом нет ничего странного. Вид уравнений состояния зависит, например, от выбора самих пере- менных состояния. Изменив набор переменных состояния, мы изменим уравне- ния состояния, хотя сама исследуемая система останется той же самой. Для ил- люстрации сказанного выполните команды: S1to2=ss2tf(S1) S2to1=tf2ss(S2) Сравните S1to2 и S2, а также S2to1 и S1. Что Вы видите? Расчёт временного отклика системы Временные характеристики — это отклики (во времени) системы на раз- личные входные воздействия. Эти воздействия могут быть типовыми или про- извольными. Среди типовых наибольшее распространение получили единичный ступенчатый сигнал и единичный импульсный сигнал. Единичный ступенчатый сигнал определяется выражением 0, при 0 1( ) 1, при 0 t t t < ⎧ = ⎨ ≥ ⎩ Единичный импульсный сигнал (единичный δ-импульс, импульс Дирака) представляет собой первую производную по времени от единичного ступенча- того сигнала. Это импульс бесконечной амплитуды и нулевой длительности, имеющий площадь, равную единице: 4 0 0, при 0 ( ) , ( ) 1 , при 0 t t d t δ δ ∞ ≠ ⎧ = = ⎨∞ = ⎩ ∫ τ τ Оба сигнала физически нереализуемы, но очень удобны для аналитиче- ского описания свойств различных объектов. Реакция системы на единичный ступенчатый сигнал называется её пере- ходной характеристикой h(t), а реакция на импульсный сигнал — импульсной характеристикой g(t). Для аналитического исследования временных характеристик используют передаточную функцию W(s). С помощью обратного преобразования Лапласа (обозначаемого символом L –1 ) находят: [ ] 1 1 W( ) ( ) ( ) W( ) s h t L s g t L s − − ⎡ ⎤ = ⎢ ⎥ ⎣ ⎦ = (6) Для расчёта временного отклика линейной системы на различные воздей- ствия в Scilab имеется универсальная функция csim(). Формат вызова: [y [,x]]=csim(u,t,sl,[x0 [,tol]]) Аргумент u определяет входное воздействие. Это может быть либо указа- тель на функцию, либо сама функция, либо специальная символьная строка: "step" — обозначает единичное ступенчатое воздействие, или "impulse" — обозначает единичное импульсное воздействие. t — вещественный вектор, задающий моменты времени для вычисления отклика. t(1) — начальное время (x0=x(t(1))). sl — линейная система (указатель на структуру syslin() ) (см. выше). y — вектор (матрица) значений выходного сигнала y=[y(t(i))], i=1,..,n x — массив значений переменных состояния x=[x(t(i))], i=1,..,n tol — вектор [atol rtol], задающий абсолютную и относительную погрешности для решателя дифференциальных уравнений. Применение функции csim() рассмотрим на следующих примерах clear //очистка всех переменных //Определим линейную систему в пространстве состояний: A=[-1, 2; -2, 0]; B=[0; 2]; C=[1, 0]; S1=syslin('c',A,B,C); //Зададим временной интервал t=0:0.02:10; //Рассчитаем переходную характеристику: h=csim("step",t,S1); 5 //Рассчитаем импульсную характеристику: g=csim("impulse",t,S1); //Построим их в графическом окне 0: scf(0);clf(0); plot2d(t',[h',g',0*t']); //Определим входной сигнал как функцию: deff("u=input1(t)","u=abs(sin(t))"); //Рассчитаем отклик на него: y1=csim(input1,t,S1); //Построим график в графическом окне 1: scf(1);clf(1); plot2d(t',[y1',input1(t'),0*t']); //Зададим входной сигнал как одиночный импульс: function u=input2(t) if t<5.0 then u=1; else u=0; end endfunction //Рассчитаем отклик на него: y2=csim(input2,t,S1); //Построим график в графическом окне 2. //Поскольку функция input2 возвращает число, а не вектор, //для построения графика входного сигнала //предварительно сформируем вектор u2(i)=input2(t(i)). n=size(t,'c'); u2=zeros(1,n); for i=1:n u2(i)=input2(t(i)); end scf(2); clf(2); plot2d(t',[y2', u2',0*t']); Копии графических окон 0, 1 и 2 приведены на рисунке 1. 6 h(t) g(t) Графическое окно 0 y1(t) Графическое окно 1 7 Графическое окно 2 y2(t) Рис. 1 Построение частотных характеристик Среди частотных характеристик наиболее известны амплитудно- частотная (АЧХ) и фазочастотная (ФЧХ) характеристики. Построенные в лога- рифмическом масштабе, они известны как логарифмические (ЛАЧХ и ЛФЧХ соответственно). Для их аналитического вычисления используется передаточ- ная функция W(s). АЧХ определяется выражением: A( ) W( ) j ω = ω , где W(j ω) получают из W(s) заменой частоты s чисто мнимой переменной jω. W(j ω) называют комплексной частотной функцией, поэтому АЧХ есть модуль комплексной частотной функции. ФЧХ есть аргумент (фазовый угол) ком- плексной частотной функции: ( ) ( ) arg W( ) j ϕ ω = ω В Scilab имеются несколько функций для исследования частотных харак- теристик линейных систем. Выше уже упоминалась функция bode(), резуль- татом которой являются графики логарифмической АЧХ и ФЧХ, построенные в одинаковом частотном масштабе. Результат вызова функции bode(S1,0.01,10) для вышеописанной системы приведён на рисунке 2. Функция выдаёт результат только в графическом виде, не инициализируя ника- ких переменных рабочей среды. Получить доступ к значениям точек АЧХ и ФЧХ можно, только обратившись к дочерним объектам графического окна. Эта процедура достаточно утомительна и неудобна. 8 Рис. 2. Диаграмма Боде Чтобы получить удобный доступ (например, программный) непосредст- венно к точкам характеристик, т.е. взять под контроль процесс построения гра- фиков, можно предложить функцию repfreq(), рассчитывающую частотный отклик линейной системы в виде последовательности значений W( ) j ω при из- менении частоты в заданном диапазоне. [[frq,] repf]=repfreq(sl,fmin,fmax [,step]) [[frq,] repf]=repfreq(sl [,frq]) [frq,repf,splitf]=repfreq(sl,fmin,fmax [,step]) [frq,repf,splitf]=repfreq(sl [,frq]) Здесь sl — указатель на структуру syslin; fmin, fmax — нижняя min f и верхняя max f граница частотного диапазо- на; frq — вещественный вектор значений частот f [Гц]; repf — вектор (комплексный) соответствующих значений частотной функции W( ) W( 2 ) j j f ω = × π ; step — шаг логарифмического приращения частоты min min min lg lg lg 2 max [10 ,10 ,10 ,..., ] f f step f step frq f + + = Если параметр step не задан, вектор frq будет вычислен с помощью встроенной функции frq=calfrq(sl,fmin,fmax). splitf — вектор индексов критических частот. Изначально splitf=[1]. Если на какой-либо частоте k f модуль частотной функции ста- 9 новится равным машинной бесконечности (т.е. не может быть вычислен), вы- числения не прерываются. Вместо этого порядковый номер этой частоты (ин- декс k) заносится в вектор индексов, увеличивая на 1 его размер: splitf=[1,k]. Размер вектора spli tf, больший единицы, говорит о том, что об ии и-либо расчётными точками наружены критические частоты. Примечание. Обнаружить критические частоты по значению вектора splitf в общем случае не удастся, поскольку нет никакой гарант , что все критические частоты совпадут с каким k f . Вероят- ность такого события ничтожно мала. Имея значения частотной функции W( ) W( 2 ) j j f ω = × π , значения АЧХ A( ω) и ФЧХ ϕ(ω) можно рассчитать по известным формулам ( ) A( ) W( ) ; ( ) arg W( ) j j ω = ω ϕ ω = ω Для этого также можно воспользоваться встроенной функцией dbphi(), которая дополнительно пересчитает АЧХ в ЛАЧХ, а фазы — в градусы: L ( ) ( ) 20lg A( ) 20lg W( ) , 180 ( ) arg W( ) j j ω = ω = ω ω π D [db,p ии; ых систем часто п ϕ ω = hi]=dbphi(repf) Здесь repf — комплексный вектор значений частотной функц db и phi — вектора значений ЛАЧХ и ФЧХ соответственно. В практике анализа линейн ользуются частотным годо- графом — частотной функцией ( ) W( ) A( ) j j e ϕ ω ω = ω , построенной на комплекс- ной плоскости. Частотный годограф представляет собой кривую, описываемую концом вектора W( ) j ω при изменении частоты, как правило, от нуля до беско- нечности. (Начало вектора закреплено в нача рдинат.) Каждое значение частоты ω определяет конкретное значение W( ) ле коо j ω и одну точку графа. Расстояние от этой точки до начала координат (т.е. длина вектора W( ) годо j ω ) даёт значение АЧХ, а угол, образованный этим вектором с положительной вещест- венной полуосью — значение ФЧХ на этой частоте. Для построения частотного годографа (также называемого АФХ или АФЧХ), можно использовать резуль- тат оп и выполните следующий пример. исанной выше функции repfreq(). Изучите clear; A=[-1,2;-2,0]; B=[0;2]; C=[1,0]; S1=syslin('c',A,B,C); //Частотный анализ 10 [f,rf]=repfreq(S1, 0.01, 9.9, 0.02); //Рассчитаем значения ЛАЧХ и ФЧХ [LA, Fi]=dbphi(rf); //Построим ЛАЧХ в графическом окне 0: scf(0); clf(0); plot2d(f,LA, style=2, logflag="ln", axesflag=1); xgrid(); xtitle("ЛАЧХ","Частота, ГЦ","20lg(A), дБ"); //Построим ЛФЧХ в графическом окне 1: scf(1); clf(1); plot2d(f,Fi, style=5, logflag="ln", axesflag=1); xgrid();xtitle("ЛФЧХ","Частота, ГЦ","Фаза, град."); //Построим АФХ (частотный годограф) в графическом окне 2: m=max([abs(real(rf)),abs(imag(rf))]); scf(2); clf(2); plot2d(real(rf),imag(rf),style=6,frameflag=3,axesflag=4, rect=[-m,-m,m,m]); xtitle("АФХ"); Копии графических окон 0, 1 и 2 приведены на рисунке 3. Графическое окно 0 11 Графическое окно 1 Графическое окно 2 Рис. 3 Графики ЛАЧХ, ЛФЧХ и АФХ Карта нулей и полюсов линейной системы Нули и полюсы линейной системы определяются по её передаточной функции: нулями называют корни полинома числителя, а полюсами — корни полинома знаменателя передаточной функции. Если передаточная функ- ция задана отношением полиномов вида i z i p 1 1 1 1 1 1 W( ) k k k k m m m m n s n s n s n s d s d s d s d − − − − + + + + = + + + 0 0 + , то через нули и полюсы её можно переписать в виде 12 1 2 1 2 ( )( )...( ) W( ) ( )( )...( ) k k m m n s z s z s z s d s p s p s p − − − = − − − , т.е. нули и полюсы с точностью до постоянного множителя определяют переда- точную функцию и соответственно, свойства системы. Нули и полюсы могут быть как вещественными, так и комплексными. Комплексные нули или полюсы всегда образуют комплексно-сопряжённые пары (с одинаковыми веществен- ными и равными по модулю, но противоположными по знаку мнимыми частя- ми). Карта нулей и полюсов представляет собой комплексную плоскость (с координатами α и jω), на которую наносятся специальными значками значения всех нулей (кружочками) и полюсов (крестиками). На их расположении осно- ваны некоторые критерии оценки устойчивости и качества системы автомати- ческого управления. Библиотечных функций для построения карты нулей и полюсов в Scilab нет. (Во всяком случае, автор этих строк их пока не обнаружил). Поэтому в ка- честве возможного варианта решения этой задачи рассмотрите следующий пример. //Определим универсальную функцию пользователя //для построения карты нулей и полюсов function KartNulPol(y) Numerator=numer(y); Denominator=denom(y); nul=roots(Numerator); pol=roots(Denominator); mx=max([abs(real(nul));abs(real(pol))]) + 0.25; my=max([abs(imag(nul));abs(imag(pol))]) + 0.25; plot2d(real(pol),imag(pol),-2,'054',rect=[-mx,-my,mx,my]); plot2d(real(nul),imag(nul),-9); //Получить маркер объектов осей текущего графич. окна a=gca(); //Задать цвет кружочков - зелёный a.children(1).children.mark_foreground=13; //Задать цвет крестиков - красный a.children(2).children.mark_foreground=5; endfunction //Определим линейную систему через передаточную функцию s=poly(0,'s'); Num=s^2+3*s+1; Den=s^3+3*s^2+3*s+10; W=syslin('c',Num,Den); //Построим карту нулей и полюсов scf(0);clf(0); KartNulPol(W); Результат выполнения приведён на рисунке 4. 13 Рис. 4. Карта нулей и полюсов Порядок выполнения работы 1. Изучить теоретические сведения, изложенные выше и выполнить приве- дённые примеры. Пользуйтесь редактором сценариев, чтобы сохранить листинги для их дальнейшей модификации. Сравнить результаты с кон- трольными графиками на рисунках 1–4. 2. Задать линейную систему Ls1 через передаточную функцию 2 2 1 0 3 2 2 1 W( ) n s n s n s 0 s d s d s d + + = + + + . Коэффициенты полиномов выбрать из табли- цы 1 по номеру бригады. Листинг команд и распечатку структуры Ls1 поместить в отчет. 3. Построить карту нулей и полюсов системы Ls1. Листинг команд и карту поместить в отчёт. 4. Преобразовать линейную систему Ls1 в Ls2, построенную в пространстве состояний. Листинг команд и распечатку структуры Ls2 поместить в от- чет. 5. Можно ли воспользоваться функцией KartNulPol(Ls2) для построения её карты нулей и полюсов? Почему? Убедиться на практике. Листинг ко- манд и результат поместить в отчет. 6. Построить переходную и импульсную характеристики системы Ls1. Оп- ределить величины их главных экстремумов. Листинг команд, графики характеристик и значения экстремумов поместить в отчёт. 7. Рассчитать отклик системы Ls1 на нетиповое воздействие. Описание входного сигнала взять из таблицы 2. Листинг команд, графики выходно- го и входного сигнала, построенные в одном графическом окне, помес- тить в отчёт. 8. Выполнить частотный анализ системы Ls2 в диапазоне частот от 0,001 до 10Гц. Листинг команд и графики ЛАЧХ, ЛФЧХ, АФХ поместить в отчёт. 14 Таблица 1. Коэффициенты полиномов передаточной функции Вариант 2 n 1 n 0 n 2 d 1 d 0 d 1. 1.0 1.10 0.100 3.0000 3.1600 1.2000 2. 1.1 1.54 0.495 2.8000 2.9200 1.2000 3. 1.2 1.08 0.096 2.3727 2.2264 0.9091 4. 1.3 1.04 0.091 2.1909 2.0264 0.9091 5. 1.4 -1.54 0.252 1.8333 1.5278 0.6944 6. 1.5 -0.90 -0.240 1.6667 1.3611 0.6944 7. 1.6 0.80 -0.224 1.3286 0.8959 0.4592 8. 1.7 1.36 0.204 1.1857 0.7673 0.4592 9. 1.8 -1.98 0.432 1.2000 0.7644 0.3556 10. 1.9 -0.76 -0.399 1.3333 0.8711 0.3556 11. 2.0 0.60 -0.360 1.2000 0.7406 0.2734 12. 2.1 1.68 0.315 1.3250 0.8281 0.2734 13. 2.2 -2.42 0.616 1.3059 0.7696 0.2076 14. 2.3 -0.46 -0.552 1.4235 0.8401 0.2076 15. 2.4 0.24 -0.480 1.3889 0.7531 0.1543 16. 2.5 2.25 0.500 1.5000 0.8086 0.1543 17. 2.6 0.26 -0.780 1.2421 0.6139 0.1108 18. 2.7 -0.27 -0.810 1.1368 0.5717 0.1108 19. 2.8 0.28 -0.840 0.8000 0.3700 0.0500 20. 2.9 3.19 0.870 0.7000 0.3500 0.0500 15 Таблица 2. Нетиповые воздействия № вариантов Воздействия 1, 6, 11, 16 2, 7, 12, 17 3, 8, 13, 18 4, 9, 14, 19 5, 10, 15, 20 Контрольные вопросы 1. Что называется передаточной функцией? 2. Что такое нули и полюсы передаточной функции? 3. Назовите и опишите типовые входные воздействия. 4. Что называется импульсной характеристикой? 5. Что называется переходной характеристикой? 6. Что называется комплексной частотной функцией? 7. Что такое АЧХ? 8. Что такое ФЧХ? 9. Что такое АФХ? 10. Могут ли матрицы 16 [ ] 1 1 0 0 2 0 2 , 1 , 1 1 , 0 0 3 4 1 ⎡ ⎤ ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ = = = − ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎣ ⎦ A B C = D определять линейную систему в пространстве состояний? Почему? 11. Как по матрице A рассчитать характеристический полином цепи? 12. Как перейти от уравнений состояния к передаточной функции? 13. Как по передаточной функции перейти к модели цепи в пространстве со- стояний? 14. Что такое диаграмма Боде? 15. Какие средства имеются в Scilab для расчета временных характеристик? 16. Какие средства имеются в Scilab для расчета частотных характеристик? 17 |