Файл-Создать-M-файлСуществуют программы двух видов: управляющие программы (Script-файлы) и процедуры (файл-функции). Все программы имеют расширение имени файла .m. Файл-функция имеет первую строку вида: Function [перечень конечных величин] = <имя процедуры> (<перечень входных величин>) Script-файлы такой строки не имеют. В файл-функциях все переменные воспринимаются как локальные, в Script-файлах все используемые переменные образуют рабочее пространство (work space), которое является единым для всех Script-файлов, вызываемых в текущем сеансе работы с системой.
Оформление М-файлов: § каждый оператор записывается в отдельной строке; в конце оператора желательно вводить символ «;» (чтобы результат действия оператора не выводился в командное окно; § можно размещать несколько операторов в одной строке, разделяя их запятой или точкой с запятой; § длинный оператор можно записывать в несколько строк, при этом предыдущая строка должна заканчиваться тремя точками (…); § строка комментариев начинается с символа %; § строки комментариев, предшествующие первому выполняемому оператору программы (должны заканчиваться пустой строкой), воспринимаются системой как описание программы. Эти строки выводятся в командное окно, если набрать в нем help <имя файла>;. § отсутствует оператор окончания текста программы; § переменные не описываются и не объявляются; их размер устанавливается при вводе значений ее элементов; § имена переменных могут содержать только буквы латинского алфавита и цифры, начинаются с буквы, число символов в имени – не больше 19-ти, строчные и прописные буквы в именах различаются системой.
Пример 2.1: Составим М-файл для вычисления значения функции: y = f1(x) = a2*ctg(x)*sqrt(sin4(x) – cos4(x). Для этого в окне редактора наберем следующий текст: function y = f1(x, a) %Процедура вычисления значения функции % y = f1(x) = a2*ctg(x)*sqrt(sin(x)^4 – cos(x)^4) %Обращение y = f1(x,a). y = (a^2)*cot(x).*sqrt(sin(x).^4 – cos(x).^4); Сохраним этот текст в файле f1.m. Введем в командное окно команду: >> y = f1(1, 0.1)
Если хотим получить вектор значений функции при различных значениях аргумента, вначале сформируем вектор аргументов: >> z = 0: 0.3: 1.8; И обратимся к процедуре следующим образом: >> v = f1 (z,1) Чтобы получить информацию о созданной процедуре, наберем >> help f1
Пример 2.2: Оформите текст примера 1.14 как файл-функцию, сопроводив ее необходимыми комментариями. Для этой процедуры нет необходимости передавать переменные из командного окна. Поэтому первую строку этой процедуры можно оформить следующим образом: function appr; Сохраните файл под именем appr.m. Вызовите программу из командного окна. Чтобы получить информацию о созданной процедуре, наберите help appr.
Пример 2.3: Измените процедуру примера 2.2 таким образом, чтобы переменные передавались из командного окна.
Некоторые стандартные файл-функции 1 [I, cnt] = quad (‘<имя функции>’, a, b) – процедура вычисления интеграла методом квадратур; здесь a, b – нижняя и верхняя границы изменения функции; I – полученное значение интеграла; cnt - число обращений к вычислению функции, представленной М-файлом с названием, указанным в <имя функции>. 2 ode23, ode45 – интегрирование обыкновенных дифференциальных уравнений; обращение к ним имеет вид; [t, y] = ode23 (‘<имя функции>’, tspan, y0, options) [t, y] = ode45 (‘<имя функции>’, tspan, y0, options) где <имя функции> - строка символов, являющихся именем М-файла, в котором вычисляется вектор-функция f(t,y), то есть правые части системы ОДУ; y0 – вектор начальных значений переменных состояния; t – массив значений аргумента, соответствующих шагам интегрирования; y – матрица матрица проинтегрированных значений переменных, в которой каждый столбец соответствует одной из переменных состояния, а строка содержит значения переменных состояния, соответствующих определенному шагу интегрирования; tspan – вектор-строка [t0 tfinal], содержащая начальное и конечное значение аргумента; options – необязательные параметры, определяющие значения допустимых относительной и абсолютной погрешности интегрирования (по умолчанию соответственно 1.0e-3 и 1.0e-6). 3 fmin – нахождение минимума функции одного аргумента; обращение: Xmin = fmin (‘<имя функции>’, X1, X2) здесь Xmin – локальный минимум функции в интервале X1<X<X2; 4 fmins - нахождение минимума функции нескольких аргументов; обращение: Xmin = fmins (‘<имя функции>’, X0), X0 – начальное значение вектора аргументов; 5 fzero - нахождение нулей (корней) функции одного аргумента; обращение; z = fzero ('<имя функции>', x0, tо1, trace) здесь x0 – начальное значение аргумента, в окрестности которого отыскивается действительный корень функции, tо1 – заданная относительная погрешность вычисления корня, trace – обозначение необходимости выводить на экран промежуточные результаты, z – значение искомого корня; 6 fplot - построение графиков функции одной переменной; отличие от процедуры plot в том, что не надо предварительно вычислять значения функции; обращение: fplot(‘<имя функции>, [<интервал>], n) где <интервал> - вектор-строка из двух чисел, задающая нижнюю и верхнюю границы изменения аргументов, <имя функции> - имя М-файла с текстом процедуры вычисления значения функции по заданному значению аргумента, n – количество частей, на которые разбит указанный интервал (по умолчанию = 25); если обратиться к этой процедуре следующим образом [x, Y] = fplot(‘<имя функции>, [<интервал>], n) то график функции не отображается на экране, а вычисляется вектор аргументов x и вектор (матрица) Y соответствующих значений указанной функции; в этом случае для построения графика следует обратиться к процедуре plot(x,Y).
Пример 2.4. Найдем число π как значение локального минимума функции y = cos(x) на отрезке [3,4]: >> Xmin = fmin(‘cos’,3,4)
Пример 2.5. Найдем минимум функции y = (x-1)^2 + 3 на отрезке [0,5]. Вначале сформируем файл-фунцию:
function y = parab (x) y = (x-1).^2 + 3 plot(x,y)
Сохраним ее в файле parab.m. Затем в командном окне наберем: >> x= -5:1:5 >>y=parab(x) >> xmin = fmin(‘parab’, 0,5)
Рассмотрим применение процедур интегрирования обыкновенных дифференциальных уравнений на следующих примерах.
Пример 2.5. Требуется решить дифференциальное уравнение y’ = -y, с начальным условием y(0)= 1. Откроем окно редактора системы и запишем в нее текст файл-функции: function dy=f1(t,y) dy=-y Сохраним файл с именем f1. В командном окне введем: >>[t, y] = ode45('f1', [0 5], 1)
Здесь мы установили интервал решения уравнения - [0 5], начальное значение переменной состояния – 1, параметр options не указывали. Система выдаст в командное окно числовые значения переменных. Чтобы просмотреть график, наберите >>plot(t, y)
Пример 2.6. Рассмотрим следующую систему:
Для решения этой системы создадим М-файл со следующим текстом: function dy = rigid(t,y) dy = zeros(3,1); % вектор-столбец dy(1) = y(2) * y(3); dy(2) = -y(1) * y(3); dy(3) = -0.51 * y(1) * y(2); Сохраним файл. Решаем систему на интервале [0 12] с вектором начальных условий [0 1 1]:
>> [t,y] = ode45('rigid',[0 12],[0 1 1]); Отразим решение на графике: >> plot(t,y)
3. Моделирование систем в пакете Simulink
Пакет SIMULINK применяется для исследования поведения динамических систем. Ввод характеристик исследуемых систем производится в диалоговом режиме, путем графической сборки схемы соединений стандартных элементарных звеньев. В результате такой сборки образуется модель исследуемой системы, которую называют S-моделью. Модель хранится в файле с расширением .mdl. Для построения S-модели используются блоки из библиотеки SIMULINK. Состав библиотеки может быть пополнен пользовательскими блоками. Библиотека блоков SIMULINK представляет собой набор визуальных объектов, используя которые можно собирать, как из кубиков, произвольную конструкцию. Для любого блока можно получать требуемое число копий и использовать каждую из них абсолютно автономно. Более того, практически для всех блоков существует возможность индивидуальной настройки: пользователь может изменить как внутренние параметры блоков (например, количество входов), так и внешнее оформление (размер, цвет, имя и т. д.). На порядок соединения блоков друг с другом также не накладывается никаких ограничений. Конечно, при связывании блоков необходимо соблюдать определенные правила, о которых будет сказано чуть позже, однако они обусловлены в основном логикой работы самой модели, а не специальными требованиями SIMULINK. Для удобства работы пользователя библиотека блоков разбита на семь разделов. Шесть из них являются базовыми и не могут изменяться пользователем (за исключением внешнего оформления): Описание состава и назначения блоков SIMULINK приведено в Приложении. Запуск пакета: Файл-Создать-Модель или соответствующая кнопка на панели инструментов. Необходимые блоки выбираются из соответствующих разделов библиотеки и с помощью мыши переносятся в окно создаваемой модели (по умолчанию это окно называется Untitled). Блоки соединяются с помощью линий, которые рисуем с помощью мыши (как от входного порта к выходному, так и наоборот). В процессе построения блок-схем все операции с объектами основаны на использовании известной технологии Drag-and-Drop (например, выделение, копирование, перестановка, отсоединение, удаление, изменение размеров, установка параметров блоков, создание соединительных линий и т.д.) Блок-схемы алгоритмов моделей составляются по содержательному описанию работы системы.
Пример 3.1. Необходимо синтезировать комбинационное цифровое устройство (КЦУ) на элементах И-НЕ на три входа, выходной сигнал которого совпадает с большинством входных сигналов. Это словесное описание условий функционирования КЦУ. Ему соответствует таблица истинности:
Соберем схему модели непосредственно по таблице истинности (рисунок 1): - в окне модели расположите три блока Constant – для входных сигналов; - преобразование входного сигнала в соответствии с заданной таблицей истинности обеспечивает блок Combinatorial Logic (раздел Nonlinear). Так как на вход этого блока может поступать только один сигнал, для объединения всех входов используйте блок Mux. Блок Combinatorial Logic имеет единственный параметр настройки — Truth table (таблица истинности), который представляет собой список возможных выходных значений автомата. При задании таблицы истинности необходимо соблюдать два основных правила: а) число строк таблицы должно быть равно 2n, где п — число элементов (размерность) входного сигнала; б) входы таблицы считаются заданными.+ - сформированный выходной сигнал выведите в блок Display.
В процессе синтеза КЦУ подразумевается необходимость минимизации аппаратных затрат на реализацию устройства. Чтобы выполнить условие задачи - синтезировать КЦУ на элементах И-НЕ на три входа, используем для моделирования совершенные нормальные формы (СНФ). Формула в СДНФ (совершенная нормальная дизъюнктивная форма): (1) После преобразований: Используя правила минимизации:
(2) Для создания модели, описывающей работу цифрового устройства на базе формулы (1), используются блоки логических операций (LogicalOperator) – рисунок 1.
Пример 3.2. Пакет MatLAB используется для исследования асинхронного двигателя в схеме асинхронного вентильного каскада (АВК). На кафедре ЭиАПУ разработаны имитационные модели структурной и электрической схем АВК[1]. На этих моделях можно выполнить серию экспериментов и получить характеристики объекта (влияние величины напряжения управления и нагрузки на характер переходных процессов, оценку соответствия установившихся значений параметров заданным и т.д.). На рисунке 2 приведена имитационная модель структурной схемы АВК, которая имеет вид последовательно соединенных звеньев. В эти звенья введены соответствующие передаточные функции с расчетными данными. Пример 3.3. Математическая модель динамических процессов в топке котельного агрегата (объект с сосредоточенными параметрами) представляется в виде системы обыкновенных дифференциальных уравнений[2]. Решение этой системы дает возможность получить закон изменения температуры газов и температуры слоя загрязнений во времени, которые могут быть использованы для регулирования процесса сжигания топлива. Рассмотрим задачу в следующей постановке: получить для различных значений расходов топлива, температур подогрева воздуха и массы загрязнений изменения температур газов и температур слоя загрязнений во времени и провести анализ влияния этих изменений на температуру поверхностей нагрева, расположенных по ходу топки. Так как часто нужно знать эти изменения при дополнительной подаче топлива или при отклонениях подогрева воздуха от номинальных значений, можно рассматривать линеаризованную математическую модель топочной камеры:
Здесь ∆t – температура топочных газов, ∆t1 – температура слоя отложений, ∆Bp – расход воздуха, ∆tB температура воздуха, τ - время, a1,…,a6 – расчетные коэффициенты модели. Структурная схема модели представлена на рисунке 3. Чтобы не загромождать диаграмму и для повышения наглядности блоки расчетов по каждому уравнению системы оформлены как подсистемы (блок Subsystem из раздела Connection). Блоки Scope позволяют в процессе моделирования наблюдать динамику изменения интересующих исследователя характеристик системы. Создаваемое с помощью блока Scope «смотровое окно» напоминает экран измерительного прибора. В подсистемах для расчета коэффициентов модели и текущих значений температур (газов и шлаков) используются блоки Fcn, Product (раздел Nonlinear), Sum (раздел Linear), Constant (раздел Sources).
4. Цифровая обработка сигналов (Signal Processing Toolbox)
Пакет Signal Processing Toolbox позволяет проектировать (рассчитывать конкретные числовые характеристики) цифровые и аналоговые фильтры по их требуемым амплитудно- и фазо-частотным характеристикам, формировать последовательности типовых временных сигналов и обрабатывать их спроектированными фильтрами. В пакет входят процедуры, осуществляющие преобразования Фурье, Гильберта, а также статистический анализ. Пакет позволяет рассчитывать корреляционные функции, спектральную плотность мощности сигнала, оценивать параметры фильтров по измеренным отсчетам входной и выходной последовательностей.
Формирование типовых процессов Процедура rectpuls обеспечивает формирование одиночного импульса прямоугольной формы. Обращение к процедуре: y = rectpuls (t, w). Здесь y – вектор значений сигнала импульса единичной амплитуды, шириной w, центрированного относительно t = 0 по заданному вектору t моментов времени. Если ширина импульса w не указана, ее значение по умолчанию принимается равным единице.
Пример 4.1: >> t = 0: 0.1: 10; >> y = 0.75*rectpuls(t-3,2) + 0.5*rectpuls(t-8, 0.4)… + 1.35*rectpuls(t-5, 0.8); >> plot(t,y), grid, set(gca, ‘FontName’, ‘Arial Cyr’, “FontSize’, 16) >> title (‘Пример применения процедуры RECTPULS’) >>xlabel (‘Время (с)’) >>ylabel (‘Выходной процесс y(t)’)
Формирование импульса треугольной формы единичной амплитуды можно осуществить при помощи процедуры tripuls, обращение к которой имеет вид: y= tripuls (t, w, s) Аргументы y, t, w имеют тот же смысл. Аргумент s (-1 < s <1) определяет наклон треугольника. Если s = 0 или не указан, треугольный импульс имеет симметричную форму.
Пример 4.2: >> t = 0: 0.1: 10; >> y = 0.75*tripuls(t-1, 0.5) + 0.5*tripuls(t-5, 0.5, -1)… + 1.35*tripuls(t-3, 0.8, 1); >> plot(t,y), grid, set(gca, ‘FontName’, ‘Arial Cyr’, “FontSize’, 16) >> title (‘Пример применения процедуры TRIPULS’) >>xlabel (‘Время (с)’) >>ylabel (‘Выходной процесс y(t)’)
Для формирования импульса, являющегося синусоидой, модулированной функцией Гаусса, используется процедура gauspuls. Если обратиться к ней по форме: y= gauspuls (t, fc, bw), то она создает вектор значений указанного сигнала с единичной амплитудой, изменяющейся с частотой fc Гц и с шириной bw полосы частот сигнала. Если два последних аргумента не указаны, они по умолчанию приобретают значения 1000 Гц и 0,5 соответственно.
Пример 4.3: >> t = 0: 0.1: 10; >> y = 0.75*gauspuls(t-3, 1, 0.5); >> plot(t,y), grid, set(gca, ‘FontName’, ‘Arial Cyr’, FontSize’, 16) >> title (‘Пример применения процедуры GAUSPULS’) >>xlabel (‘Время (с)’) >>ylabel (‘Выходной процесс y(t)’)
Процедуры sin(x), cos(x) формируют колебания, состоящие из конечного числа гармонических составляющих.
Пример 4.4: >> t = 0: 0.01: 50; >> y = 0.7*sin(pi * t/5); >> plot(t,y), grid, set(gca, ‘FontName’, ‘Arial Cyr’, “FontSize’, 16) >> title (‘Гармонические колебания y = 0.7*sin(pi * t/5) >>xlabel (‘Время (с)’) >>ylabel (‘Выходной процесс y(t)’)
Процесс, являющийся последовательностью прямоугольных импульсов с периодом 2π для заданной в векторе t последовательности отсчетов времени, генерируется при помощи процедуры square. Обращение: y = square (t, duty), где аргумент duty определяет длительность положительной полуволны в процентах от периода волны.
Пример 4.5: >> y = 0.7*square(pi*t/5, 40); >> plot(t,y), grid, set(gca, ‘FontName’, ‘Arial Cyr’, “FontSize’, 16) >> title (‘Прямоугольные волны y = 0.7*square(pi*t/5, 40); >>xlabel (‘Время (с)’) >>ylabel (‘Выходной процесс y(t)’)
Генерирование пилообразных и треугольных колебаний выполняется процедурой sawtooth. Обращение: y= sawtooth (t, width) Параметр width определяет часть периода, в которой сигнал увеличивается.
Пример 4.6: >> y = 0.7* sawtooth (pi*t/5, 0.5); >> plot(t,y), grid, set(gca, ‘FontName’, ‘Arial Cyr’, “FontSize’, 16) >> title (‘Прямоугольные волны y = 0.7* sawtooth (pi*t/5, 0.5); >>xlabel (‘Время (с)’) >>ylabel (‘Выходной процесс y(t)’)
Процедура pulstran позволяет формировать колебания, являющиеся последовательностью прямоугольных, либо треугольных, либо гауссовых импульсов. Обращение y = pulstran (t, d, ‘func’, p1, p2, …) Здесь d определяет вектор значений тех моментов времени, где должны быть центры соответствующих импульсов; параметр func определяет форму импульсов; параметры p1, p2,… определяют необходимые параметры импульса в соответствии с формой обращения к процедуре, определяющей этот импульс.
Пример 4.7: >.> t= 0: 0.01: 50; >> d = [0: 50/5: 50]; >> y = 0.7* pulstran (t, d, ‘tripuls’, 5); >> plot(t,y), grid, set(gca, ‘FontName’, ‘Arial Cyr’, FontSize’, 16) >> title (‘y = 0.7* pulstran (t, d, ‘tripuls’, 5’); >>xlabel (‘Время (с)’) >>ylabel (‘Выходной процесс y(t)’)
Аналогично сформируйте последовательность прямоугольных и гауссовых импульсов.
Общие средства фильтрации Для контроля и графического представления передаточной функции любого линейного динамического звена удобно использовать процедуру freqs. Обращение к ней: h = freqs (b, a, w) Процедура создает вектор h комплексных значений частотной характеристики W(jω) по передаточной функции W(s) звена, заданной векторами коэффициентов ее числителя – b и знаменателя - a, а также по заданному вектору w частоты ω. Если аргумент w не указан, процедура автоматически выбирает 200 отсчетов частоты, для которых вычисляется частотная характеристика. Если не указана выходная величина, процедура выводит в текущее графическое окно два графика - АЧХ и ФЧХ.
Пример 4.8: Передаточная функция имеет вид W (s) = a/(s2 + 2*ξ*ω0*s + ω02 ) Выведем графики АЧХ и ФЧХ. >> T0 =1; dz =0.05; >> om0 = 2*pi/T0; A = 1; >>a1(1) = 1; a1(2) = 2*dz*om0; a1(3) = om0^2; b1(1)=A; >>freqs (b1,a1)
Чтобы получить дискретную частотную характеристику по дискретной передаточной функции, заданной векторами значений ее числителя b и знаменателя a, используется процедура y = freqz (b, a)
Спектральный анализ процессов Функции fft (Fast Fourier Transformation) и ifft (Invers Fast Fourier Transformation) осуществляют преобразования заданного вектора, соответствующие дискретному прямому и обратному преобразованию Фурье. Обращение: y = fft (x, n); x = ifft (y, n), где n - число элементов заданного вектора (а также и выходного вектора).
Пример 4.9: >> t = 0: 0.001:2; >> x = sin(2*pi*5*t) + cos(2*pi*12*t); %Формируем входной процесс >> plot (t, x); grid >> set (gca, ‘FontName’, ‘Arial Cyr’, ‘Fontsize’, 16); >> title (‘Входной процесс’); >> xlabel (‘Время’); >> ylabel (‘x(t)’) >> y = fft (x); >> a = abs (y); >> plot (a); grid >> set (gca, ‘FontName’, ‘Arial Cyr’, ‘Fontsize’, 16); >> title (‘Модуль Фурье-изображения’); >> xlabel (‘Номер элемента вектора’); >> ylabel (‘abs(F(x(t)’)
Выполним обратное преобразование.
Пример 4.10: >> x = ifft (y); >> plot (t, x); grid >> set (gca, ‘FontName’, ‘Arial Cyr’, ‘Fontsize’, 16); >> title (‘Обратное Фурье-преобразование); >> xlabel (‘Номер элемента вектора’); >> ylabel (‘x(t)’)
Для определения массива Фурье-изображения с целью построения графика Фурье-изображения в частотной области используется процедура fftshift. Пример 4.11: >> f1 = -500: 0.5:500; >> v = fffshift (y); >> a = abs (v); >> plot (f1(970:1030), a(970:1030)); grid >> set (gca, ‘FontName’, ‘Arial Cyr’, ‘Fontsize’, 16); >> title (‘Модуль Фурье-изображения’); >> xlabel (‘Частота (Гц)’); >> ylabel (‘abs(F(x(t)’)
По графику спектра невозможно определить амплитуды гармоник. Для этого нужно вектор Фурье-изображения разделить на количество его элементов, амплитуды распределяются между положительными и отрицательными частотами поровну, поэтому они вдвое меньше истинной амплитуды соответствующей гармоники.
Пример 4.12: >> N = length (y); >> a = abs (v)/N; >> plot (f1(970:1030), a(970:1030)); grid >> set (gca, ‘FontName’, ‘Arial Cyr’, ‘Fontsize’, 16); >> title (‘Модуль Фурье-изображения’); >> xlabel (‘Частота (Гц)’); >> ylabel (‘abs(F(x(t)/N’)
Фильтрация осуществляется процедурой filter следующим образом: y = filter (b, a, x), где x – заданный вектор входного сигнала; y - вектор значений выходного сигнала фильтра, получаемого в результате фильтрации; b - вектор коэффициентов числителя дискретной передаточной функции фильтра; a - вектор коэффициентов знаменателя этой функции. Эта функция обеспечивает формирование вектора y по заданным b, a, x, в соответствии с конечно-разностным уравнением фильтра с дискретной передаточной функцией вида рациональной дроби:
y(k) = b(1)*x(k) + b(2)*x(k-1) + b(nb+1)*x(k-nb)- - a(2)*y(k-1) – a(3)*y(k-3) -… - a(na+1)*y(k-nb), где b=[b(1),…,b(nb+1)], a=[1,a(2),…,a(nb+1)].
Пример 4.13: Пусть требуется получить достаточно верную информацию о некотором полезном сигнале (Yp) синусоидальной формы с периодом T1 = 1 и амплитудой A1 = 0.75, к которому добавился шум первичного преобразователя в виде более высокочастотной синусоиды с периодом T2 = 0.2 и амплитудой A2 = 5, а в результате измерения еще добавился белый гауссовый шум измерителя с интенсивностью Aш = 5. Сформируем этот сигнал как вектор его значений в дискретные моменты времени с дискретом Ts = 0.001. В результате получим такой измеренный сигнал x(t):
>> Ts = 0.001; >> t =0: Ts:20; >> A1 =0.75; T1 = 1; >> T2 = 0.2; A2 = 10; eps = pi/4; >> Ash = 5; >> x = A1.*sin92*pi*t./T1 + A2.*sin(2*pi.*t./T2+eps) + … Ash*randn(1, length(t)); >> plot (t(10002:end),x(10002:end), grid, set(gca,‘FontName’, … ‘ArialCyr’, ‘FontSize’, 16), >> title(‘Входной процесс’); >> xlabel (‘Время (с)’); >> ylabel (‘X(t)’)
Требуется так обработать измеренные данные x(t), чтобы восстановить по ним полезный сигнал как можно точнее. Период собственных колебаний фильтра должен быть равен периоду колебаний полезного сигнала Tф = T1. Чтобы после прохождения через фильтр амплитуда восстановленного сигнала совпадала с амплитудой полезного сигнала, нужно входной сигнал фильтра домножить на постоянную величину 2ξω02. Сформируем этот фильтр, пропустим через него сформированный процесс, для сравнения на одном графике отобразим и восстанавливаемый процесс:
Пример 4.14: >> T1 = 1; Tf = T1; dz = 0.05; >> om0=2*pi/Tf; A=1; oms = om0*Ts; >> a(1)=1+2*dz*oms+oms^2; >> a(2)=-2*(1 + dz*oms); >> a(3)= 1; >> b(1)=A*Ts*Ts*(2*dz*om0^2); >> Yp=A1*sin(2*pi*t/T1); >> y=filter(b,a,x); >> plot(t(1002:end), y(t(1002:end),t(1002:end), Yp(1002:end), grid,… set (gca,’,‘FontName’, ‘ArialCyr’, ‘FontSize’, 16), >> title(‘Процесс на выходе фильтра'); >> xlabel (‘Время (с)’); >> ylabel (‘Y(t)’)
Чтобы избежать фазовых искажений полезного сигнала при его восстановлении, можно воспользоваться процедурой двойной фильтрации – filtfilt.
Пример 4.15: >> y = filtfilt (b,a,x) >> plot(t, y, t, Yp, grid, set (gca,’,‘FontName’, ‘ArialCyr’, ‘FontSize’, 16) >> title(‘Применение процедуры filtfilt’); >> xlabel (‘Время (с)’); >> ylabel (‘Y(t)’)
Статистический анализ процессов Определение спектральной плотности выполняется процедурой psd. Обращение: [S,f] = psd (x, nfft, Fmax), где x – вектор заданных значений процесса, nfft - число элементов вектора, которые обрабатываются процедурой fft, Fmax = 1/Ts - значение частоты дискредитации сигнала, S - вектор значений СП сигнала, f – вектор значений частот, которым соответствуют найденные значения СП. В общем случае длина двух последних векторов равна nfft/2.
Пример 4.16: >> [C,f] = psd (y1, dovg, Fmax); >> subplot >> stem (f(1:200), C(1:200)), grid, >> set (gca,’,‘FontName’, ‘ArialCyr’, ‘FontSize’, 16), >> title (‘Спектральная плотность'); >> xlabel (‘Частота (Гц)’);
Если эту же процедуру вызвать без указания выходных величин, то результатом ее выполнения станет выведение графика зависимости СП от частоты.
Пример 4.17: >> psd (y1, dovg, Fmax)
Здесь значения СП откладываются в логарифмическом масштабе в децибелах.
Группа функций xcorr вычисляет оценку взаимной корреляционной функции (ВЛФ) двух последовательностей x и y. Обращение: c= xcorr (x, y) вычисляет и выдает вектор c длины 2N-1 значений ВКФ векторов x и y длины N. Также позволяет вычислить автокорреляционную функцию (АКФ) последовательности, заданной в векторе x.
Пример 4.18: >> R = xcorr (y1); >> tau = -10+Ts: Ts: 10; >> lt = length (tau); >> s1r = round (length (R) /2 – lt/2; >> s2r = round (length (R) /2 + lt/2 -1; >> plot (tau, R (s1r:s2r), grid >> set (gca,’,‘FontName’, ‘ArialCyr’, ‘FontSize’, 16), >> title (‘АКФ случайного процесса'); >> xlabel (‘Запаздывание (с)’);
Мы не рассматривали здесь вопросы проектирования фильтров в MatLab, а также интерактивную оболочку SPTool, содержащую средства поиска и просмотра сигналов (Signal Browser), проектировщик фильтров (Filter Designer), средство просмотра характеристик фильтров (Filter Viewer), средство просмотра спектра (Spectrum Viewer).
5. Исследование линейных стационарных систем (пакет Control Toolbox)
Пакет позволяет работать с процедурами анализа и синтеза линейных стационарных систем автоматического управления.
Ввод и преобразование моделей Ввод модели линейной стационарной системы (русскоязычное сокращение ЛСС, здесь же LTI - линейные, инвариантные во времени системы) в среду пакета CONTROL возможен в трех формах – в форме матриц пространства состояния, в виде коэффициентов числителей и знаменателей передаточных функций и в форме задания нулей, полюсов и коэффициента передачи системы. Некоторые процедуры создания LTI-моделей: ss - создает модель пространства состояния по заданным матрицам уравнений состояния системы; dss - создает аналогичную модель по описанию пространства состояния более общего вида, когда уравнения переменных состояния не разрешены относительно производных; tf - создает модель по заданным передаточным функциям системы zpk - создает модель по заданным нулям, полюсам и коэффициентам передачи системы; filt - создает модель по дискретным передаточным функциям, записанным в форме полиномов от z—1. set - присваивает значения некоторым другим полям –объектов (названиям входов и выходов, названиям системы и т.д.). Указанные процедуры позволяют создавать как непрерывные, так и дискретные модели. Также они применяются для преобразования моделей из одной из указанных форм в другую.
Пример 5.1: >> v = tf ([1 4], [1 2 100])
В MatLab предусмотрена возможность «набирать» программно «схему» САУ путем предварительного ввода моделей звеньев и последующего «соединения» этих звеньев в единую структуру. Однако более простым и удобным средством создания систем из отдельных блоков является интерактивная система Simulink, рассмотренная в п.3. Чтобы получить отдельные характеристики (матрицы и векторы, описывающие пространство состояния, коэффициенты числителя и знаменателя передаточной функции и т.д.) полученной модели используются процедуры: tfdata - для получения векторов числителя и знаменателя передаточной функции системы; ssdata - для получения значений матриц уравнений пространства состояния; zpkdata - для получения векторов значений полюсов и нулей системы. get - дает возможность получить полную характеристику модели.
Пакет предоставляет широкий набор процедур, осуществляющих анализ САУ с самых различных точек зрения и, прежде всего, определение откликов системы на внешние воздействия как во временной, так и в частотной области: impulse - нахождение отклика на единичное импульсное входное воздействие; step - нахождение реакции системы на единичный скачок входного воздействия; initial - определение собственного движения системы при произвольных начальных условиях; lsim - определение реакции системы на входное воздействие произвольной формы, задаваемое в виде вектора его значений во времени. Кроме того: процедуры представления реакции системы на внешние гармонические воздействия в частотной области (bode, nyquist, nichols, sigma); процедуры, вычисляющие отдельные характеристики и графически показывающие расположение полюсов и нулей системы (pole, zpkdata,gram, damp, pzmap, rlocus). Если в командном окне набрать ltview, то можно вызвать окно обозревателя LTI–объектов и в интерактивном режиме строить все графики, причем для нескольких систем одновременно.
Некоторые процедуры, используемые для синтеза САУ: lqr, lqry - осуществляют проектирование линейно-квадратичного оптимального регулятора для систем непрерывного времени; lqrd, dlqr - для проектирования дискретного оптимального линейно- квадратичного регулятора; kalman - выполняет расчет фильтра Калмана для непрерывных или дискретных САУ.
Список литературы 1. Гультяев А, Имитационное моделирование в среде Windows. – С-Пб.: “КОРОНАпринт”, 1999. 2. Лазарев Ю. MatLAB 5.x. – Киев: «Ирина», BHV, 2000.
Содержание
ПРИЛОЖЕНИЕ БИБЛИОТЕКА МОДУЛЕЙ (БЛОКОВ) SIMULINK Библиотека блоков SIMULINK представляет собой набор визуальных объектов, используя которые можно собирать, как из кубиков, произвольную конструкцию. Для любого блока можно получать требуемое число копий и использовать каждую из них абсолютно автономно. Более того, практически для всех блоков существует возможность индивидуальной настройки: пользователь может изменить как внутренние параметры блоков (например, количество входов), так и внешнее оформление (размер, цвет, имя и т. д.). На порядок соединения блоков друг с другом также не накладывается никаких ограничений. Конечно, при связывании блоков необходимо соблюдать определенные правила, о которых будет сказано чуть позже, однако они обусловлены в основном логикой работы самой модели, а не специальными требованиями SIMULINK. Для удобства работы пользователя библиотека блоков разбита на семь разделов. Шесть из них являются базовыми и не могут изменяться пользователем (за исключением внешнего оформления): • Sources (Источники), • Sinks (Получатели), • Discrete (Дискретные элементы), • Linear (Линейные элементы), • Nonlinear (Нелинейные элементы), • Connections (Соединения). Седьмой раздел — Blocksets &Toolboxes (Наборы блоков и инструменты) — содержит блоки, относящиеся к компонентам MATLAB, включенным пользователем в рабочую конфигурацию пакета. При этом для каждой компоненты создается свой подраздел библиотеки. При минимальной рабочей конфигурации в разделе Blocksets &Toolboxes имеется только один подраздел - SIMULINK Extras (Дополнение к SIMULINK). Этот подраздел, в свою очередь, разбит на шесть частей, три из которых являются дополнением к одноименным основным разделам библиотеки, а три других имеют самостоятельное значение. Это наборы блоков Transformations (блоки пересчета координат и шкал температуры), Flip-Flops (блоки, соответствующие основным типам триггеров): Linearization (содержит единственный блок, реализующий функцию линейной аппроксимации). Каждый блок, входящий в библиотеку SIMULINK, имеет по крайней мере один параметр настройки. Задавая требуемое значение параметра (или выбирая его из предлагаемого меню), пользователь имеет возможность скорректировать функцию, реализуемую данным блоком. Чтобы открыть окно настройки параметров, нужно дважды щелкнуть на изображении блока. Однако возможность изменять значения параметров появляется только после того, как блок будет помещен в окно блок-диаграммы. Окна настройки параметров всех библиотечных блоков имеют идентичную структуру и содержат краткую характеристику блока, поля ввода (или выбора) значений параметров блока и 4 кнопки: • Apply — применить; • Revert — вернуть предыдущее значение параметров; • Help — вызов файла помощи в формате html; • Close — закрыть окно настроек. Измененные значения параметров вступают в силу после «нажатия» кнопки Apply. Чтобы запустить модель на исполнение с новыми параметрами, закрывать окно настроек не обязательно. Вернемся к основным разделам библиотеки SIMULINK. Как уже было сказано, их шесть. Раздел Sources (Источники) Блоки, входящие в этот раздел, предназначены для формирования сигналов, обеспечивающих управление работой модели в целом или отдельных ее частей - блоки-источники имеют по одному выходу и не имеют входов. Итак, в качестве источников сигналов (входных величин) могут использоваться следующие блоки: Constant — формирует постоянную величину (скаляр, вектор или матрицу); Signal Generator — создает непрерывный сигнал произвольной формы; Step — генерирует единичный дискретный сигнал с заданными параметрами; Ramp — создает линейно возрастающий (убывающий) сигнал; Sine Wave — генератор гармоничес
|
