close

Вход

Забыли?

вход по аккаунту

?

1927 borodenko v.a matematicheskie zadachi energetiki i kompyuternoe modelirovanie v.a.borodenko

код для вставкиСкачать
Министерство образования и науки
Республики Казахстан
Павлодарский государственный университет
им. С. Торайгырова
Кафедра Электрические станции и автоматизация энергосистем
МЕТОДИЧЕСКИЕ УКАЗАНИЯ
к лабораторному практикуму по дисциплине
“Математические задачи энергетики
и компьютерное моделирование”
для студентов специальностей 21.01.40,21 03.40
всех форм обучения
П авлодар
УДК 51+004.94]:621.311(07)
ББК 22.1+32.973J:31.2*7
М34
Рекомендовано ученым советом института энергетики
и автоматизации ПГУ им. С. Торайгырова
Рецензенты:
канд. техн. наук, доцент, Амбарников Г.А.
Бороденко В.А.
Математические задачи энергетики и компьютерное моделиро­
вание. Лабораторный практикум. - Павлодар, Изд-во ПГУ, 2004. 27 с.
В методических указаниях к лабораторному практикуму приво­
дятся рекомендации по выполнению лабораторных работ к дисцип­
лине «Математические задачи энергетики и компьютерное моделиро­
вание», исходная схема моделируемой системы и варианты заданий,
изложены последовательность выполнения работ и порядок оформ­
ления отчета. Моделирование линейных систем базируется на
предварительном изучении теории автоматического управления.
В качестве среды для компьютерного моделирования линейных
систем используется матричная лаборатория MalLAB (версия 6)
фирмы The MathWorks, Inc.
С. Торайг
атыидагы
академик С.БеЛсембаев
атындагы гыльчии
К 1Т А П Х А Н Д С Ы
О
С Павлодарский государственный университет им. С. Торайгырова,
2004
ш , at
й ъ ч
Содержание
1
2
3
4
Задание
Визуально-ориентированное моделирование систем
Моделирование систем во временной области
Моделирование систем в частотной области
Коррекция характеристик моделируемых систем
Литература
4
5
10
16
22
26
Задание
1 Структурная схема системы автоматического регулирования
(САР) скорости двигателя постоянного тока
Рисунок 1
Двигатель (блок 6) описывается дифференциальным уравнением
r r i d 2y
L
ь
Т — ^-+ 2Т ,с— + у = к,х.
3 dt
* dt 7 ^
2 Исходные данные
Таблица 1
ki
кг
кз
*4
Г,
Тг
Т3
Ъ
Ва эиант(равен номеру студента по порядку в списке группы)
1,11 2,12 3,13 4, 14 5, 15 6, 16 7. 17 8, 18 9, 19 10,20
1,00 1.10 1,20 1,30 1,40 1,50 1,60 1.70 1.80 1,90
12,00 11,00 10,00 9,00 8,00 7,00 6,00 5,00 4,00 3,00
равен номеру по порядку первой буквы фамилии в алфавите
блок 4 соер инен с точкой а (вариант 1-10) или Ь (вариант 11-20)
3,00 4,00 5,00 6,00 7,00 8,00 9,00 10,00 11,00 12,00
1,00 2,00 3,00 4,00 5,00 6,00 7,00 8,00 9,00 10,00
0,10 0,20 0,30 0,40 0,50 0,60 0,70 0,80 0,90 1,00
0,50 0,80 1,10 1,40 1.70 1,50 1,20 1,00 0,80 0,60
0,90 0,80 0,80 0,60 0,60 0,40 0,40 0,20 0,20 0,15
Примечание - коэффициент кос обратной связи (блок 5) равен
единице. Отчет к каждой лабораторной работе выполняется отдельно
в текстовом редакторе Word в соответствии с требованиями ГОСТ
2.105-95, печатается на принтере и сдается преподавателю. Каждая
работа начинается с названия и описания цели работы.
4
1 Визуально-ориентированное моделирование систем
1.1 Цель работы
Используя заданную структурную схему и данные своего вари­
анта, создать в среде Simulink матричной лаборатории MATLAB гра­
фическую модель САР, получить численные параметры моделей сис­
темы в виде переменных состояния, передаточной матрицы, описания
нулями и полюсами, оценить устойчивость системы.
1.2 Указания к работе
Запустить программу MATLAB с помощью ярлыка на рабочем
столе или открыв файл \MATLAB6pl\BIN\WIN32\matlab.exe. Оста­
вить для работы только командное окно (Command Window), закрыв
остальные окна. Через меню File-New-Model открыть окно конструи­
рования модели. Кнопкой Simulink или Library Browser на панели ин­
струментов открыть окно библиотеки элементов.
Перемещая мышью блоки Передаточная функция (Transfer Fen раздел Continues), Усилитель Gain и Сумматор Sum (раздел Math) из
окна элементов в окно конструирования, соединяя их входы и выхо­
ды, составить структурную схему системы (см. пункт 1.3). Командами
Flip и Rotate в контекстном меню Format изменить, если нужно, ори­
ентацию блока. Кроме того, следует поменять стандартные названия
блоков на указанные в задании.
Открывая двойным щелчком левой кнопки или правой кнопкой
мыши свойства элемента (Block parameters), задать численные вели­
чины. Коэффициенты полиномов записывают в виде векторов в квад­
ратных скобках, последовательно, начиная с коэффициента при стар­
шей степени s, разделяя пробелами (саму переменную s не пишут). У
инвертирующего входа сумматора заменяют плюс на минус. Если в
изображении блока видны не все коэффициенты, нужно увеличить
размер блока (растянуть прямоугольник).
Выбрав из набора блоков Control System Toolbox элементы Вход
(Input Point), связать первый из них с входом r(t), второй с входом f(t),
элементы Выход (Output Point) связать аналогично с выходами y(t),
e(t), дав соответствующие названия. Сохранить созданную модель на
диске в папке своей группы, задав имя в формате RPxxYYzz (для спе­
циальности 21.03) и PSxxYYzz (для специальности 21.01), где х х - две
цифры номера группы, YY - две цифры номера варианта, zz - две
цифры номера первой буквы фамилии. Расширение mdl будет при­
своено файлу автоматически.
5
Через меню Tools-Linear Analysis (Инструменты-Линейный ана­
лиз) следует вызвать окно графического представления результатов
моделирования (характеристик) LTI Viewer. Далее командой меню
Simulink-Get Linearized Model (Заимствовать модель) передают чис­
ленные значения математической модели системы графическому ана­
лизатору. Автоматически будут построены четыре переходные харак­
теристики - h ef(t) и h er (t) в окнах первой строки, hyf(t) и hyr(t)b
окнах второй строки. Наконец, через меню File-Export-Export to Work­
space следует передать данные в рабочее пространство (базу
данных) Workspace программы MATLAB.
Если теперь в командном окне ввести после приглашения »
имя модели в том виде, как оно было показано в графическом окне (с
дополнительным индексом), то будет распечатано основное матема­
тическое описание системы матрицами А, В, С, D в пространстве со­
стояний. При этом строки матриц соответствуют выходам, а столбцы
- входам системы или отдельных блоков.
Рекомендуется далее избавиться от дополнительного индекса,
переименовав модель, например, командой (имя использовать свое)
»w w O 10234=ww010234_1;
где слева требующееся имя системы, а справа - присвоенное ей гра­
фическим анализатором (регистр букв учитывается). Затем откройте
окно рабочего пространства через меню View-Workspace и удалите
все объекты, кроме своей системы ww010234. Сохраните данные на
диск с помощью кнопки Save Workspace As..., задав нужное имя фай­
ла (расширение mat программа добавит автоматически) и выбрав пап­
ку своей группы или папку Work в каталоге MatLAB.
Теперь имеется модель системы в виде структурной схемы
(файл .mdl) и набор присущих ей параметров (файл .mat). В заключе­
ние работы получим модели системы в виде пространства состояний
(State Space), передаточных функций (Transfer Function), нулей, полю­
сов и коэффициентов усиления (Zero-Pole-Gain) командами
»ss(wwO10234)
% в скобках использовать имя своей системы
»tf(w w 010234)
»zpk(ww010234)
Коэффициент усиления системы в установившемся режиме
(матрица коэффициентов для многомерной системы) может быть по­
лучен командой
»dcgain(wwO10234)
Многомерная система устойчива, если все её собственные зна­
чения (корни характеристического полинома) имеют отрицательную
действительную часть. Для вычисления собственных значений систе­
6
мы можно использовать функцию eig(mui системы) или последова­
тельность команд
% вычисляем характеристический полином системы
»ро1у(имя системы)
»roots(ans)
% вычисляем корни полинома
Осталось перенести графическую модель из среды Simulink ко­
мандой Edit-Сору to ClipBoard и полученные описания из командного
окна (команды Сору-Paste) через буфер в отчет по лабораторной рабо­
те (окно Word), сопроводив их необходимыми пояснениями.
1.3 Пример оформления отчета
Структурная схема системы wwO10234, построенная в среде
Simulink
к ю
Описание системы в пространстве состояний
»ss(ww 010234)
а=
sl/W5
sl/W5
sl/W3
S1/W2
S1AV5
-10
1
0
0
F(t)
S1AV5
S1/W5
sl/W3
S1/W2
sl/W5
-10
0
0
-20
R(t)
1
0
0
0
0
0
0
1
7
sl/W3
1
0
-1
-20
sl/W2
0
0
1
-1
E(t)
Y(t)
sl/W5 sl/W5
0
-20
0
20
sl/W3
0
0
sl/W2
0
0
d=
E(t)
Y(t)
F(t)
0
0
R(t)
1
0
Continuous-time model.
Описание системы передаточными функциями
»tf(ww010234)
Transfer function from input "F(t)" to output...
-20 sA2 - 40 s - 420
sA4 + 12 sA3 + 51 sA2 + 230 s + 230
20 sA2 + 40 s + 420
Y(t):------------------------------sA4 + 12 sA3 + 51 sA2 + 230 s + 230
Transfer function from input "R(t)" to output...
sA4 + 12 sA3 + 51 sA2 + 230 s + 210
E(t):------ -----------------------sA4 + 12 sA3 + 51 sA2 + 230 s + 230
20
sA4 + 12 sA3 + 51 sA2 + 230 s + 230
Описание системы нулями, полюсами и коэффициентом переда» zpk(ww010234)
Zero/pole/gain from input HF(t)” to output...
-20 (^2 +2s+21)
m . ------------------------------
(s+8.841) (s+1.258) (syv2 + 1.901s + 20.68)
20 (sA2 + 2s + 21)
Y(t): -------------------------------------(s+8.841) (s+1.258) (s/42 + 1.901s + 20.68)
Zero/pole/gain from input "R(t)” to output...
(s+8.873) (s+1.127) (sA2 + 2s + 21)
E(t): -------------------------------------(s+8.841) (s+1.258) (s*2 + 1.901s + 20.68)
20
Y(t): -------------------------------------(s+8.841) (s+1.258) (sA2 + 1.901s + 20.68)
Коэффициент усиления системы в установившемся режиме
» dcgain(ww010234)
ansт
-1.8261 0.9130
1.8261 0.0870
Проверка устойчивости системы
» eig(ww010234)
ans “
-8.8412
-1.2379
-0.9505 + 4.4472i
-0.9505 - 4.4472i
Все собственные значения системы (корни характеристического
полинома) имеют отрицательную действительную часть, поэтому сис­
тема устойчива.
9
2 Моделирование систем во временной области
2.1 Цель работы
Используя созданную модель, изучить принципы моделирова­
ния поведения системы и получить её основные характеристики во
временной области - переходную функцию (Step Response) с прямыми
показателями качества, импульсную функцию (Impulse Response), ре­
акцию на начальные условия (Response to Initial Conditions) и разло­
жение на простые дроби (Residue).
2.2 Указания к работе
Откройте командное окно программы MATLAB, загрузите ра­
бочее пространство с данными своей модели (файл .mat).
Внимание! Если система с исходными значениями параметров
оказалась неустойчивой, сейчас её следует сделать устойчивой, под­
бирая значение коэффициента \ц в структурной схеме и повторяя дей­
ствия из первого раздела до получения нужной переходной характе­
ристики новой системы - у устойчивой системы свободная состав­
ляющая со временем затухает. При этом исключите в контекстном
меню графика (Systems) предыдущие варианты системы, чтобы они не
мешали наблюдать измененные характеристики. Сохраните на диске
новую модель и рабочее пространство, укажите новое значение 1с«и
причину его изменения в отчете.
Переходная характеристика (Step Response) строится по умол­
чанию после вызова графического анализатора командой
LTIVIEW(HMa модели). Автоматически появятся четыре окна (по чис­
лу входов и выходов системы) с графиками h(t) переходной функции
по соответствующему каналу.
Проведем анализ только для главной передаточной функции от входа задания R(t) до выхода управляемой величины Y(t), что соот­
ветствует элементу W ^ s) во второй строке и втором столбце матри­
цы передаточных функций системы. Для этого, вызвав щелчком пра­
вой кнопки мыши на любом графике контекстное меню, выберем I/O
Selector, а в его окне щелкнем нижний правый маркер и закроем окно.
Определим параметры оценок качества, получаемых из пере­
ходной характеристики. Для этого, вызвав контекстное меню графика,
выберите Properties (Свойства) и в окне редаюирования свойств вкладку Characteristics (Характеристики). Задайте показ максимума
(Peak Response), времени регулирования (Settling Time), времени на­
растания (Rise Time) и установившегося значения (Steady State или
DC Gain). Зону А для оценки времени регулирования следует устано-
10
вить равной 5 %, интервал оценки времени нарастания можно оста­
вить прежним - от 10 до 90 %.
Устанавливая поочередно указатель мыши на маркеры, записать
значения времени нарастания tu, времени регулирования tp, коэффици­
ента усиления в установившемся режиме h(<x>), максимального значе­
ния характеристики h— времени достижения максимума
и пере­
регулирования о (Overshoot) в процентах. Скопировать график снача­
ла в буфер, выбрав опции File-Print to Figure, а в окне Figure опции
Edit-Сору Figure, затем из буфера в отчет (окно Word - вставить).
Использовав в контекстном меню графика команды Plot ТуреImpulse, получим импульсную характеристику системы, обозначаемую
обычно w(t) или g(t).
Для получения переходной и импульсной временных характери­
стик можно также использовать команды step(wM системы), grid и
impulse(HM* системы), grid, при этом сразу выводится окно Figure.
Вычисление реакции на начальные условия производится MatLAB только для систем, описанных в пространстве состояний, по за­
данному вектору начальных условий хО. Учтем при расчете, что полу­
чению импульсной характеристики соответствует начальное значение
переменной состояния на входе системы (обычно это х«), равное еди­
нице. Поэтому для функции передачи Wy,(s), которая равна элементу
W(2,2) матрицы передаточных функций системы, необходимо выпол­
нить команды
»х0=[0 0 0 1];
% формируем вектор начальных значений
»initial(ww010234(2,2),хО), grid
В результате получают график реакции на начальные условия,
совпадающий по форме с импульсной характеристикой (численно
значения могут отличаться из-за особенностей структурной схемы).
Реакцию системы в аналитическом виде (оригинал) по ее изо­
бражению Лапласа можно найти с помощью функции residue, задав
полиномы Ъ, а числителя и знаменателя рациональной дроби изобра­
жения
»[r,pjc]=residue(b,a)
Векторы коэффициентов Ь, а можно ввести самостоятельно, а
можно получить из существующей модели, например, командой
»[b,a]=tfdata(ww010234(2,2),V)
Функция residue выводит значения вычетов (коэффициентов
числителей), полюсов и свободной части разложения на простые дро­
бя. Вычеты (элементы вектора г) печатаются в том же порядке, в ка­
ком выводятся полюса (вектор р) - первый сверху вычет соответству­
ет первому полюсу, т.е. корню характеристического уравнения и т.п.
11
Пользуясь таблицей соответствий оригиналов и их изображений по
Лапласу, для действительного полюса Xс вычетом а можно записать
оригинал в виде экспоненты
^ + ± а * е х р ( ± Я * /) >
± А»
а для комплексной пары полюсов и вычетов (здесь а, Р - действитель­
ная и мнимая части вычета, X, со - действительная и мнимая части по­
люса)
± а ± jP
±Л ± ja>
можно записать
2 * ехр(± X * t) *(± а * cos(tw * t) ± р * sin(o> * f))
(плюс перед Р записывается при несовпадении знака мнимой части у
соответствующих друг другу вычета и полюса).
Имея аналитическую запись оригинала, можно рассчитать и по­
строить временную характеристику, например, организовав цикл типа
FOR ... END в интервале времени от 0 до
для заданного количест­
ва точек (время
лучше всего взять из графика переходной характе­
ристики). Длинную формулу можно продолжить на другой строке,
прервав знаком троеточия. В формуле для комплексных сопряженных
корней умножение записывается с точкой впереди, так как произво­
дится умножение массивов, а не матриц. Точка с запятой в конце опе­
ратора отменяет его немедленное выполнение и вывод результата.
2.3 Пример оформления отчета
Передаточная функция системы Wyi(s) равна
» wyr=tf(wO10234(2,2))
Transfer function from input "R(t)" to output "Y(t)":
20
sA4 + 12 sA3 + 51 ^ 2 + 230 s + 230
Переходная характеристика hy^t) с прямыми оценками качества
12
»hiview(wyr)
Из графика время нарастания равно 1,74 с, время регулирования
2,36 с, время достижения максимума 4,08 с, значение максимума
0,0868, перерегулирование 0 %, установившееся значение 0,087.
Импульсная характеристика guW имеет вид
»impulse(wyr), grid
tapin Response
РгавЯОО
13
Реакция многомерной системы WW010234 на начальные усло­
вия х1(0)=0; х2(0)=0; х3(0)=0; х4(0)=1 по выходу Y(t) имеет вид
» х 0 = [0 0 0 1];
»initial(ww010234(2,2),x0), grid
Response to HM Concttons
Полученная характеристика по форме совпадает с импульсной.
Найдем оригинал передаточной функции W ^ s) с помощью раз­
ложения на простые дроби. Коэффициенты полиномов числителя и
знаменателя функции-изображения
»[b,a]=tfdata(w010234(2,2), V)
b=
0
0
0
0 20.0000
а=
1.0000 12.0000 51.0000 230.0000 230.0000
Коэффициенты числителей и корни знаменателей разложения на
простые дроби
» [r,p,k]=residue(b,a)
14
-0.0321
-0.0503 + 0.0239i
-0.0503 - 0.0239i
0.1327
P=
-8.8412
-0.9505 + 4.4472i
-0.9505 - 4.4472i
-1.2579
к=
□
Расчет значений оригинала в интервале времени 0 - tycr, с с ша­
гом 0.01 с и построение графика (система имеет два действительных
полюса и пару комплексных сопряженных)
»t=0:0.01:4.5;
»y=-0.0321*exp(-8.8412*tyH).1327*exp(-1.2579*t)+...
2*exp(-0.9505*t).*(-0.0503*cos(4.4472*t)-0.0239*sin(4.4472*t));
» plot(t,y), grid;
» ЬЙеСОригинал передаточной функции*);
Оригинал передаточной функции
Оригиналом передаточной функции является импульсная функ­
ция, поэтому полученная характеристика по виду и числовым значе­
ниям совпадает с найденной ранее импульсной.
15
3 Моделирование систем в частотной области
3.1 Цель работы
Используя созданную модель, изучить принципы моделирова­
ния системы в частотной области и получить её основные характери­
стики - амплитудно-фазовую характеристику (АФЧХ) или годограф
Найквиста (Nyquist Diagram), логарифмические характеристики (ЛЧХ)
или диаграммы Боде (Bode Diagram), запасы устойчивости системы,
амплитудную (АЧХ), фазовую (ФЧХ), действительную (ВЧХ) и мни­
мую (МЧХ) характеристики.
3.2 Указания к работе
Откройте командное окно программы MATLAB, загрузите ра­
бочее пространство с данными своей модели (файл .mat). Будем ана­
лизировать только коэффициент передачи от входа R(t) к выходу Y(t),
что соответствует элементу передаточной матрицы Wm(s)
»wyr=tf(ww010234(2,2))
Если, вызвав графический анализатор командой ltiview(wyr),
выбрать в контекстном меню графиков тип графика Plot Type-Nyquist,
будет построена АФЧХ. Напомним, что MATLAB строит обе кривых
- как основную в диапазоне частот со от 0 до плюс бесконечности, так
и симметричную ей в диапазоне частот от 0 до минус бесконечности,
основная кривая определяется по направлению стрелки на ней и зна­
чениям частот.
Если необходим анализ устойчивости, то следует задавать для
анализа передаточную функцию разомкнутой системы, а на графике
оценивать положение особой точки с координатами (-1, jO) - АФЧХ
неустойчивой в замкнутом состоянии системы охватывает эту точку.
Если устойчивость не анализируется, рекомендуется изменить мас­
штаб графика командой Zoom так, чтобы вписать в него кривую пол­
ностью независимо от положения точки (-1, jO).
Логарифмические амплитудная (ЛАЧХ) и фазовая (ЛФЧХ) час­
тотные характеристики строятся по команде Plot Type-Bode (с сеткой
Grid на графике, если необходимо). Отметим, что ЛАЧХ (Magnitude)
вычисляется по формуле L(w)=201gA(w) и измеряется в децибелах
(dB), ЛФЧХ (Phase) откладывается в градусах, а частота (Frequency)
по оси абсцисс откладывается в рад/с и в логарифмическом масштабе.
По команде Оатр(имя системы) выводятся величины собствен­
ных значений (Eigenvalue) или полюсов системы, соответствующие
им значения коэффициента демпфирования или затухания %
(Damping) и частот сопряжения, а по команде Dcgain(mui системы) -
16
значение коэффициента усиления на нулевой частоте (добротность),
что бывает полезно при самостоятельном построении JIA4X.
Функция Evalfir позволяет вычислить комплексный отклик
(Frequency Response) системы в виде действительной RealO и мнимой
ImagO составляющих на единственный гармонический сигнал задан­
ной частоты, например
» fr=0+j *0.5;
% задаем комплексную частоту
» evalfr(w(2,2),fr)
% получаем действительную и мнимую часть
реакции
ans =
0.0723 - 0.0378i,
а амплитуду и фазу комплексного отклика можно получить с помо­
щью функций AbsO и AngleO
Графики нетиповых функций обычно строят либо с помощью
команды Plot, либо используя команды semilogxO и semilogyO - с ло­
гарифмическим масштабом только по оси х или только по оси у,
loglogO - с логарифмическим масштабом по обеим осям графика.
Диапазон частот можно задавать в виде вектора f=0:0.1:10, в ли­
нейном f=linspace(0,10,100) или логарифмическом f=(-2,2,100) мас­
штабе, причем во всех трех случаях было задано 100 точек. В первых
двух командах - для диапазона частот от 0 до 10 рад/с, в последней
команде для диапазона частот от 10'2 до 10+2 рад/с.
Запасы устойчивости системы по амплитуде и фазе обычно оп­
ределяются для разомкнутой системы с помощью логарифмического
критерия Найквиста, для чего используется команда margin
»margm(cncTeMa), grid
при этом запас по амплитуде Gm указывается в децибелах, запас по
фазе Pm - в градусах. Для обеих величин фиксируется частота (рад/с)
и отмечается положение на графике. Значение Inf говорит о том, что
по данному параметру система имеет бесконечный запас устойчиво­
сти, отрицательные или нулевые значения - об отсутствии запасов.
Передаточную функцию разомкнутой системы можно получить
одним из следующих способов:
а) вычесть соответственно из коэффициентов полинома знаменателя
передаточной функции замкнутой системы значения коэффициентов
полинома числителя, оставив полином числителя неизменным;
б) в блоке 5 модели системы изменить значение кос на ноль, т.е. разомкнуть систему, после чего уже извецшымтяам путем ^раздел 1)передать параметры модели в рабочее пространство MATLAB. Р
атында*ы ПМУ-дИ
(академик С.Бейсембаев
атындагы гылыми
I / п » г- v д и Л Г К1
1 K I I м* lA A rlM V ^D I
17
Первый способ может применяться только в случае охвата сис­
темы главной единичной отрицательной обратной связью, второй бо­
лее универсален.
3.3 Пример оформления отчета
Передаточная функция системы Wy^s) равна
» wyr=tf(wO10234(2,2))
Transfer function from input ”R(t)" to output "Y(t)":
20
sA4 + 12 sA3 + 51 sA2 + 230 s + 230
Диаграмма Боде
»bode(w yr)
Собственные значения (полюса), коэффициенты демпфирования
и частоты собственных колебаний системы
» damp(wyr)
Eigenvalue
-1.26е+000
-9.50е-001 + 4.45e+000i
-9.50е-001 - 4.45e+000i
-8.84е+000
Damping Freq. (rad/s)
1.00е+000
1.26e+000
2.09e-001
4.55e+000
2.09e-001
4.55e+000
1.00e+000
8.84e+000
18
Вычисляем реакцию системы Y на единичное гармоническое
воздействие в линейном диапазоне частот от 0 до 10 рад/с (200 точек),
переводя действительную частоту в комплексную
» f=linspace(0,10,200);
» for i=l :200
y(i)=evalfr(wyr,0+j*f(i));
end
% задаем вектор частот
% организуем цикл
Совмещаем вещественную (синий цвет - ‘blue’) и мнимую
(красный цвет - ‘red’) характеристики на одном графике
»
»
»
»
»
»
hold on;
% фиксируем графическое окно
plot(f,real(y),,b');
% выводим график в окно
plot(f;imag(y),V);
grid;
titleCReal and Imaginary Characteristics');
hold off
ншмand HiM jm vy Characteristics
01
0.0B
0.04
0 02
-0.02
-0.04
■0.08
0
1
2
3
4
5
6
7
6
9
Строим амплитудную и фазовую характеристики системы в
двух окнах одновременно
»
»
»
»
»
»
subplot(2,l,l);
plot(f,abs(y)), grid;
title(,A4X');
ylabel('Амплитуда');
subplot(2,l,2);
plot(f,angle(y)*57.7), grid;
% задаем первое окно
% рисуем АЧХ
% делаем подписи
% задаем второе окно
% рисуем ФЧХ с переводом в градусы
19
» titleC<D4X');
» у1аЬе1('Фаза, град')
АЧС
Амплитудно-фазная частотная характеристика (подписи на гра­
фике и оси добавлены в режиме Edit окна Figure)
» plot(real(y),imag(y)), grid;
» title ('АФЧХ')
АФЧХ
Разомкнем систему, вычтя коэффициенты полинома числителя
из коэффициентов полинома знаменателя передаточной функции
»[Ъ , al^fcU ^w yr/v’);
»a=a-b;
»wol=tf(b,a)
после чего можем оценить устойчивость замкнутой системы по обеим
формам критерия Найквиста.
20
Годограф Найквиста (АФЧХ разомкнутой системы)
»nyquist(wol), grid
АФЧХ разомкнутой системы не охватывает точку с координа­
тами (-1, jO), поэтому после замыкания система будет устойчивой.
Запасы устойчивости системы по амплитуде и фазе (логарифми­
ческий критерий Найквиста)
»margm(wol), grid
ВойМдп
O n - 26022 сВ («С4.3771 radfcec). P m -hf
После замыкания запас системы по амплитуде равен 26.022 дБ
(на частоте 4,377 рад/с), запас по фазе равен бесконечности (180°).
21
4 Коррекция характеристик моделируемых систем
4.1 Цель работы
Используя созданную модель, изучить принципы моделирова­
ния системы нулями-полюсами на комплексной плоскости, метод
корневого годографа (Root Locus), методы коррекции параметров и
синтеза одномерных (SISO) и многомерных (MEMO) систем с задан­
ными свойствами.
4.2 Указания к работе
Откройте командное окно программы MATLAB, загрузите ра­
бочее пространство с данными своей модели (файл .mat). Будем ана­
лизировать либо всю систему, либо только коэффициент передачи от
входа R(t) к выходу Y(t), что соответствует элементу передаточной
матрицы W22(s) и одномерному представлению системы
»wyr=tf(wwO10234(2,2))
Для систем, описанных в пространстве состояний, обычно оце­
нивают наблюдаемость и управляемость. Для этого командами obsvO
и ctrbO формируются матрицы наблюдаемости и управляемости, а за­
тем с помощью функции гапк() вычисляется ранг матрицы. Напомним,
что система полностью управляема (наблюдаема), если ранг матрицы
управляемости (наблюдаемости) равен порядку системы.
Вызвав окно графического анализатора командой ltiview(wyr),
выберите в контекстном меню графиков тип графика Plot ТуреPole/Zero, в результате этого будет выведена карта размещения нулей
Zero (кружки) и полюсов Pole (крестики) системы на комплексной
плоскости. Аналогичные действия выполняет функция ргтар(имя
системы). Отдельно вычисляет полюса, т.е. корни полинома знамена­
теля передаточной функции системы, функция pole(), нули системы,
т.е. корни полинома числителя, - функция zeroOФункция RlocusO производит расчет и вывод в графическое ок­
но корневого годографа (Root Locus). Корневым годографом называ­
ется совокупность траекторий движения на комплексной плоскости
корней полинома знаменателя R(s) = D(s) + k*N(s) передаточной
функции замкнутой системы при изменении положительного вещест­
венного коэффициента к от 0 до бесконечности. Здесь D(s) - знамена­
тель передаточной функции разомкнутой системы, N(s) - её числи­
тель. Начало всех ветвей годографа находится в полюсах разомкнутой
системы, часть ветвей заканчивается в нулях разомкнутой системы
(если таковые имеются), остальные уходят в бесконечность.
22
Для каждой точки годографа можно получить с помощью мы­
ши: значение коэффициента передачи Gain, значение полюса замкну­
той системы Pole, значения коэффициента демпфирования Damping и
частоты собственных колебаний Frequency (рад/с), значение перерегу­
лирования Overshoot (%). У неустойчивой системы действительная
часть хотя бы одного полюса положительна, а значение коэффициента
демпфирования этого корня отрицательно.
Функция гкоо1(имя системы) открывает окно SISO Design Tool.
Она позволяет спроектировать методом корневого годографа для од­
номерной системы вводимый дополнительно в схему регулятор
(Compensator), наблюдая результирующее изменение переходной ха­
рактеристики, годографа Найквиста или диаграммы Боде (меню
Tools). Регулятор С при этом может располагаться в цепи как прямой
(Forward), так и обратной (FeedBack), как отрицательной, так и поло­
жительной связи (кнопка +/-). Место установки регулятора выбирает­
ся с помощью кнопки FS, получаемая в результате передаточная
функция регулятора выводится в окне. Цель коррекции характеристик
системы выбирается студентом совместно с преподавателем.
В процессе моделирования могут быть исключены, компенсиро­
ваны или добавлены действительные и комплексные полюса или ну­
ли, заданы характеристики звена Prefilter (фильтра F) на входе задания
и измерительного преобразователя Sensor (датчика Н) в цепи обрат­
ной связи. Исходный объект (Plant) сохраняется неизменным в виде
блока G. В меню Tools может быть выбран вид исследуемой характе­
ристики (с помощью опции Loop Responses), либо построена специ­
ально для целей исследования структурная схема в среде Simulink.
4.3 Пример оформления отчета
Передаточная функция Wy^s) равна
» wyi=tf(w010234(2,2))
Transfer function from input "R(t)n to output "Y(t)n:
20
sA4 + 12 sA3 + 51 sA2 + 230 s + 230
Проверка управляемости исходной системы
» ctrb(ww010234)
23
ans -
1.0000 0 -10.0000
0 90.0000 1.0000 -800.0000 -12.0000
О
0 1.0000
0 -10.0000
0
90.0000 1.0000
О О О
1.0000
0
-2.0000 -20.0000 -17.0000
О
1.0000 0
-1.0000 -20.0000-19.0000 220.0000 59.0000
» rank(ans)
ans =
4
Ранг матрицы управляемости равен порядку системы 4, поэтому
система полностью управляема
Проверка наблюдаемости исходной системы
» obsv(ww010234)
ans =
1.0e+003 *
0
-0.0200
0
0
0.0200
0
0
0
-0.0200
0
0
0
0.0200
0
0
0
0.2000 0.2000 -0.0200
0
-0.2000 -0.2000 0.0200
0
-1.8000 -2.0000 0.2200 -0.0200
1.8000 2.0000 -0.2200 0.0200
» rank(ans)
ans =
4
Ранг матрицы наблюдаемости равен порядку системы 4, поэто­
му система полностью наблюдаема
Нули (отсутствуют) и полюса моделируемой одномерной сис­
темы W22(S)
» zero(wyr)
ans =
Empty matrix: 0-by-l
» pole(wyr)
ans =
24
-8.8412
-0.9505 + 4.4472i
-0.9505 - 4.4472i
-1.2579
Корневой годограф системы W22(s)
»rlocus(wyr)
Root Locus
Окно проектирования методом корневого годографа системы с
заданными параметрами
25
Л итература
1 Электрические системы. Математические задачи электроэнер­
гетики: Под ред. А. Веникова. Учебник, 3-е изд. - М.: Высшая школа,
1991.-238 с.
2 Лукас В. А. Теория автоматического управления: Учеб. для
вузов. - 2-е изд., перераб. и доп. - М.: Недра, 1990. - 416 с.
3 Лазарев Ю. Ф. MatLAB 5.x - К., Издательская группа BHV,
2000.-384 с.
4 Гультяев А. Визуальное моделирование в среде MATLAB:
учебный курс. - СПб: Питер, 2000. - 432 с.
5 Дьяконов В., Круглов В. MATLAB. Анализ, идентификация и
моделирование систем. - СПб.: Питер, 2002. - 448 с.
6 Васильков Ю. В., Василькова Н. Н. Компьютерные техноло­
гии вычислений в математическом моделировании: Учеб. пособие. М.: Финансы и статистика, 2001. - 256 с.
26
Документ
Категория
Без категории
Просмотров
1
Размер файла
569 Кб
Теги
modelirovanie, 1927, borodenko, energetic, zadachi, matematicheskih, kompyuternye
1/--страниц
Пожаловаться на содержимое документа