close

Вход

Забыли?

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

?

Mironovskiy 0BBF6341AC

код для вставкиСкачать
МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РОССИЙСКОЙ ФЕДЕРАЦИИ
Федеральное государственное автономное
образовательное учреждение высшего образования
САНКТ-ПЕТЕРБУРГСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ
АЭРОКОСМИЧЕСКОГО ПРИБОРОСТРОЕНИЯ
Л. А. Мироновский
МОДЕЛИРОВАНИЕ ЛИНЕЙНЫХ СИСТЕМ
в MATLAB и SIMULINK
Лабораторный практикум
Санкт-Петербург
2017
УДК 004.4
ББК 32.965.4
М64
Рецензенты:
кандидат технических наук, доцент Г. С. Бритов;
доктор технических наук, профессор Л. А. Осипов
Утверждено
редакционно-издательским советом университета
в качестве лабораторного практикума
Мироновский, Л. А.
М64 Моделирование линейных систем в MATLAB и SIMULINK:
лабораторный практикум / Л. А. Мироновский. – СПб.: ГУАП,
2017. – 74 с.
Содержатся теоретический материал и практические примеры
к выполнению лабораторных работ по компьютерному моделированию с использованием программных пакетов MATLAB и SIMULINK.
Предназначено для проведения лабораторных работ по курсу «Моделирование» со студентами дневного обучения по направлению
09.03.01 «Информатика и вычислительная техника».
УДК 004.4
ББК 32.965.4
© Мироновский Л. А., 2017
© Санкт-Петербургский государственный
университет аэрокосмического
приборостроения, 2017
ПРЕДИСЛОВИЕ
Целью выполнения лабораторных работ является закрепление
теоретического материала по курсу лекций и получение практических навыков компьютерного моделирования.
На каждую лабораторную работу отводится два занятия. На первом проводится коллоквиум, во время которого проверяется качество написания отчета, теоретические знания студентов, а также
знание пакета программ, в котором выполняется работа. На втором
занятии работа выполняется на компьютере. Работа считается выполненной, если результаты компьютерных экспериментов хорошо
совпадают с теоретическими. Отметка, которая ставится в журнал,
учитывает качество отчета, ответы на коллоквиуме и результаты
практического выполнения работы.
Отчет оформляется каждым студентом индивидуально. Отчет
должен содержать постановку задачи, вывод необходимых формул,
расчеты, графики, программы или схемы моделирования. При выполнении расчетов в построении графиков следует придерживаться
следующих правил.
1. Для построения графика функции y = f(x) или y = f(t) сначала
составляется таблица. Полезно также найти особые точки функций, производную в начале координат, период и т. д., что облегчает
построение графика.
2. При вычислениях следует удерживать три значащие цифры.
3. При построении графиков масштабы по осям следует выбирать
в соответствии с диапазоном изменения переменных.
4. Число точек для построения графика зависит от вида функции. Например, для функции y = e–t/τ достаточно взять 5–6 точек на
интервале 0 ≤ t ≤ 3τ (τ – постоянная времени). Для периодических
кривых надо брать несколько точек на полпериода, особенно тщательно рассчитав первый период. На графике должно поместиться
1–3 периода. Выбор шага при компьютерном моделировании производится из аналогичных соображений.
3
Методические указания содержат описания семи лабораторных работ, которые выполняются в программной среде MATLAB
и SIMULINK.
В конце приведен необходимый справочный материал – математические формулы, таблицы преобразований Лапласа и z-преобразований, значения функций e–t, sint, cost и некоторые команды MATLAB.
4
Лабораторная работа № 1
МОДЕЛИ ЛИНЕЙНЫХ БЛОКОВ
Цель работы: ознакомиться с характеристиками и методикой компьютерного моделирования блоков с дробно-рациональными передаточными функциями с помощью пакетов MATLAB
и SIMULINK.
1. Теоретические сведения
1.1. Математическое описание линейных блоков
Рассмотрим линейный блок Q, преобразующий входной сигнал
x(t) в выходной сигнал y(t) (рис. 1).
Напомним, что блок называется линейным, если он удовлетворяет принципу суперпозиции, который нестрого формулируется
следующим образом: следствие от суммы причин равно сумме следствий от каждой из причин, взятых в отдельности.
Это означает, что если y1 и y2 - реакции блока на входные сигналы x1 и x2, то реакция блока на входной сигнал x = c1x1 + c2x2
должна быть равна y = c1y1 + c2y2 при любых числах c1, c2 и любых
сигналах x1, x2 .
Например, линейность интегратора, описываемого уравнением y = a∫xdt, следует из известного правила «интеграл суммы равен
сумме интегралов»:
y=
a ∫ (c1x1 + c2 x2 )dt =
c1a ∫ x1dt +c2a ∫ x2dt =
c1y1 + c2 y2 .
Для квадратора, описываемого уравнением y = x2, принцип суперпозиции не выполняется, так как y = (x1 + x2)2 ≠ x12 + x22 , поэтому квадратор является нелинейным блоком.
Удобным средством для моделирования структурных схем в системе MATLAB является SIMULINK. Он содержит библиотеки различных блоков, из которых в рабочем окне с помощью мышки строится структурная схема модели.
x
Q
y
Рис. 1. Исследуемая система
5
Рассмотрим некоторые из типовых вычислительных блоков.
Масштабный усилитель (обозначается треугольником с надписью gain) реализует математическую зависимость y = ax, где а – коэффициент усиления. Если а = –1 то масштабный усилитель превращается в инвертор – блок, изменяющий знак входного сигнала.
Инвертор реализует зависимость у = – х.
Сумматор (обозначается кружком, треугольником или прямоугольником с надписью sum) реализует зависимость y = ±x1 ± x2 ± x3,
знаки суммирования и число слагаемых можно изменять, щелкнув
по блоку мышкой. Оба эти блока находятся в библиотеке (группе
блоков) Math Operations.
Простейший динамический блок – интегратор (изображается
1
треугольником или прямоугольником с обозначением
внутри
s
и надписью integrator). Он находится в библиотеке Continuos и реализует зависимость
t
y(t=
) y0 + ∫ x(t) dt.
0
Для установки начального условия у0 надо войти внутрь блока,
щелкнув по нему мышкой. Например, если у0 = 2 и на вход интегратора поступает сигнал х = 1, то выходной сигнал определяется формулой y = 2 + t. Это линейно изменяющееся напряжение, осциллограмма которого имеет вид наклонной прямой.
Соединяя блоки определенным образом между собой, получают
схемы для реализации различных функций времени, решения дифференциальных уравнений, моделирования технических объектов.
В частности, для получения экспоненциальных функций можно
использовать интегратор, охваченный положительной или отрицательной обратной связью (рис. 2, а, б).
Чтобы найти вид выходного сигнала первой из этих схем, выпишем уравнение для ее выходного сигнала
dy
=
y ay
=
; y(0) y0 , где y =
.
dt
Это дифференциальное уравнение первого порядка относительно
переменной y. Решим его методом разделения переменных:
dy
=
a dt; ln y =
at + C;
y
y=
y0 eat .
Аналогичные выкладки для второй схемы приводят к выражению y = y0 e–at. Графики этих сигналов показаны на рис. 2, в.
6
y0
а)
y0
б)
y
y
y
a
–a
в)
y
y
а
y0
б
t
Рис. 2. Получение экспоненциальных функций
В библиотеке Continous имеется также блок для моделирования
произвольной передаточной функции (transfer function). Он обозначается прямоугольником с надписью Transfer Fcn и обозначением
1
внутри. Он, в частности, позволяет моделировать линейные
s +1
блоки, которые можно описать дифференциальным уравнением
первого порядка
a1y (t) + a0 y(t) = b1x (t) + b0 x(t). (1)
При исследовании и моделировании бортовых систем такими
уравнениями пользуются для описания динамики автопилотов,
сервоприводов и измерительных датчиков. Постоянные коэффициенты a0, a1, b0, b1 зависят от технических параметров соответствующих устройств.
Линейному блоку вида (1) соответствует передаточная функция
первого порядка.
Определение 1. Передаточной функцией Q(p) линейного блока
называется отношение изображений по Лапласу его выходного и
входного сигналов
Q( p) =
Y ( p)
X( p)
(2)
при нулевых начальных условиях.
7
Напомним, что изображение по Лапласу F(p) функции f(t) задается формулой
∞
F ( p) = ∫ f (t)e- pt dt. (3)
0
Оно определено для функций f(t), равных нулю при t < 0.
Пример 1. Функцией единичного скачка (единичной функцией)
называется функция
0 ïðè t < 0,
x=
(t) 1=
(t) 
1 ïðè t ≥ 0.
Ее изображение по Лапласу имеет вид
∞
∞
1
1
- pt
X( p) =
- e- pt =
.
∫ e dt =
p
p
0
0
Производная по времени от функции единичного скачка называd
ется дельта-функцией: d(t) = 1(t). Она равна нулю везде, кроме
dt
точки t = 0, в которой она принимает бесконечно большое значение.
Ее инженерная интерпретация – очень короткий импульс единичной площади. Операции дифференцирования по времени соответствует умножение на оператор р в области изображений по Лапласу,
1
поэтому для изображения функции d(t) получаем p = 1.
p
В пакете MATLAB для получения единичной функции имеется
команда stepfun. В SIMULINK эта функция получается с помощью
блока step, находящегося в библиотеке Sources (источники).
Установим вид передаточной функции блоков первого порядка.
Для этого применим к обеим частям уравнения (1) преобразование
Лапласа (начальные условия считаем нулевыми)
a1 pY (p) + a0 Y (p) = b1 pX(p) + b0 X(p).
b1 p + b0
Y (p)
= Q=
(p)
.
X(p)
a1 p + a0
Здесь X(p) и Y(p) - изображения сигналов x(t) и y(t) по Лапласу,
Q(p) – передаточная функция (ПФ) блока Q.
Замечание 1. В зарубежной литературе (и в SIMULINK) в качестве аргумента передаточной функции используется буква s. Мы
употребляем букву р, придерживаясь традиций отечественной
учебной литературы.
Отсюда получаем
8
Замечание 2. В определении ПФ не оговаривается вид входного
сигнала. Дело в том, что для линейных блоков отношение Y(p)/X(p)
не зависит от вида x(t) (для отношения y(t)/x(t) это, разумеется, неверно). Поэтому при нахождении ПФ можно использовать любой
входной сигнал, не равный тождественно нулю.
Пример 2. Найдем передаточную функцию интегратора.
Полагая x(t) = 1(t), получим y(t) = t, t≥0. Изображения входного
и выходного сигналов, соответственно, равны 1/p и 1/p2, следовательно, ПФ интегратора имеет вид Q(p) = 1 / p.
Пример 3. Найдем передаточную функцию усилителя с коэффициентом усиления k. Полагая x(t) = 1(t), получим y(t) = k · 1(t), т. е.
k 1
ПФ Q=
( p) =
:
k.
p p
1.2. Прохождение сигналов через линейные звенья
Определение ПФ с помощью формулы (2) остается справедливым
для линейных стационарных звеньев произвольного порядка, а также для схем, составленных из таких звеньев.
Удобство использования передаточной функции состоит в том,
что она позволяет определять реакцию y(t) звена Q на любой конкретный входной сигнал x(t). Для этого находят изображение X(p)
входного сигнала (по таблице преобразований Лапласа) и умножают его на передаточную функцию Q(p), получая тем самым изображение выходного сигнала Y(p). Затем, используя таблицу, выполняют обратный переход от найденного изображения к оригиналу y(t).
Пример 4. Найдем реакцию сервопривода с передаточной функцией Q(p) = k/(Tp + 1) (апериодическое звено) на управляющее воздействие x(t) = 1(t). Коэффициент усиления сервопривода k = 1, постоянная времени T = 0,1 c.
1
В соответствии с примером 1 X( p) = , следовательно,
p
1
Y ( p) =
. Представим Y(p) в виде суммы простейших дробей
p(0,1 p + 1)
1
0,1
1
1
Y ( p) =
=
.
p 0,1 p + 1 p
p + 10
Обращаясь к таблице преобразований Лапласа, находим соответствующий оригинал y(t) = 1–e–10t. Графики сигналов x(t) и y(t)
приведены на рис. 3. Величина τ = 3T приблизительно характеризует время отработки входного сигнала сервоприводом.
9
x
1
y
t
Т
2Т
3Т
Рис. 3. Переходная функция сервопривода
Пример 5. Найдем реакцию изодромного звена с передаточной
kp
функцией Q( p) =
, на входной сигнал X = cos(ωt), если k = 2,
Tp + 1
–
1
T = 1 c, w = 1 c .
p
Изображение входного сигнала имеет вид X( p) = 2
, поэтому
p +1
2p
2 p2
p
Y ( p) =
⋅ 2
=
.
p + 1 p + 1 ( p + 1)( p2 + 1)
Выполняем разложение на простейшие дроби
2 p2
A
Bp
C
.
=
+ 2
+ 2
( p + 1)( p + 1) p + 1 p + 1 p + 1
2
Постоянные A, B, C находим, приводя выражение в правой части к общему знаменателю и приравнивая коэффициенты
при одинаковых степенях p в числителях правой и левой части:
A = B = 1; С = –1. Возвращаясь к оригиналам, получаем y(t) = e–t +
+ cos t – sint.
Ввиду того, что реакция звена на импульсное или ступенчатое воздействие исчерпывающим образом характеризует линейное
звено, в теории управления для них используются специальные
термины.
Определение 2. Импульсной весовой характеристикой (весовой функцией) q(t) называется реакция звена на входной сигнал
x(t) = d(t). Импульсной переходной характеристикой (переходной функцией) p(t) называется реакция звена на входной сигнал
x(t) = 1(t).
10
d
Эти функции связаны соотношением q(t) =
p(t), которое слеdt
d
дует из равенства d(t) = 1(t). Отсюда вытекает, что при компьюdt
терном моделировании весовую характеристику блока с ПФ Q(p),
т. е. реакцию на x(t) = d(t), можно получить, подавая единичный
скачок на блок с ПФ pQ(p). Другой способ – использовать сигнал
с входа выходного интегратора схемы.
2. Задание по работе и содержание отчета
В первой части лабораторной работы исследуются схемы моделирования, показанные на рис. 4, а, б, в. Числовые значения парамета)
5
б)
5
x
–1
x
–1
1
y0
y0
–k
–k
y
в)
y
5
x
–1
z
y0
y
–k
с
Рис. 4. Исследуемые схемы моделирования
11
ра k приведены в таблице вариантов, начальные условия и значе1
ние с рассчитываются по формулам y0 =
-10 + 15k, c =
.
2 - 3k
Отчет по первой части работы должен содержать:
1. Схемы рис. 4 с заданными численными значениями y0, k, c, нарисованные применительно к SIMULINK (с указанием блока входного сигнала и осциллографов).
2. Вывод формул для выходных сигналов каждой из схем. Таблицы и графики этих сигналов, координаты точки экстремума функции z(t) и тангенс угла ее наклона в начале координат.
Во второй части лабораторной работы исследуется линейное звено первого порядка, которое описывается дифференциальным уравнением y + ay =
bx. Параметры а, b, берутся из таблицы вариантов.
В отчете требуется для заданных значений а и b найти теоретически
выходной сигнал блока y(t), если входной сигнал имеет вид:
а) x(t) = 1(t), б) x(t) = e – t, в) x(t) = sin t, г) x(t) = d(t).
Отчет по второй части должен содержать:
1. Дифференциальное уравнение и передаточную функцию заданного блока.
2. Вывод формул y(t) для всех вариантов входного сигнала.
3. Таблицы и графики теоретических значений y(t) для всех случаев.
4. Схемы моделирования в пакете SIMULINK.
3. Порядок выполнения работы
Первая часть
1. С помощью пакета SIMULINK собрать схему рис. 4, а, установить коэффициенты усиления и начальные условия. Наблюдать
выходные сигналы схемы с помощью осциллографа Scope. Зарегистрировать с помощью осциллографа момент пересечения графиков
x(t) и y(t), сравнить его с теоретическим.
2. Собрать схему рис. 4, б. Для этого в предыдущей схеме (а) отключить входной сигнал и охватить каждый интегратор обратной
связью. Наблюдать выходные сигналы интеграторов с помощью осциллографа, сравнить их с теоретическими.
3. Собрать схему рис. 4, в. Для этого взять сумматор, и усилитель
с коэффициентом с и соединить их входы с выходами интеграторов
из предыдущей схемы (б). Наблюдать осциллограмму сигнала z(t),
12
сравнить теоретические и экспериментальные координаты точки
экстремума.
Вторая часть
1. Построить схему моделирования, содержащую генератор
входного сигнала и блок с заданной передаточной функцией.
2. Выбрать параметры (численный метод, время моделирования)
и выполнить моделирование. Наблюдать результат моделирования
в виде графика.
3. Получить переходные функции звеньев второго порядка с передаточными функциями 1/(p2 + b), p/(p2 + b), 1/(p2 + ap + b),
1/(p2 – ap + b).
4. Выполнить те же пункты в пакете MATLAB, используя команды step, impulse, lsim, plot.
4. Контрольные вопросы
1. Какие блоки, объекты и системы называются линейными? Являются ли линейными следующие блоки: усилитель; интегратор;
апериодическое звено; усилитель, включенный последовательно
с диодом; пассивная интегрирующая RC-цепочка?
2. Имеется интегратор с нулевыми начальными условиями. Нарисовать вид выходного сигнала, если на вход интегратора подаются:
а) х = 1(t), б) x = sin t, в) x = соs t,
г) x = et, д) x = е– t
е) х = t, ж) х = t2.
3. Выполнить задание п. 2 для случая двух одинаковых последовательно включенных интеграторов.
4. Как подать выходной сигнал интегратора на вход осциллографа? Как изменить длительность развертки осциллографа? Как наблюдать одновременно два сигнала?
5. Найти передаточную функцию усилителя, интегратора, апериодического звена, двух последовательно соединенных интеграторов. Написать соответствующие дифференциальные уравнения.
6. Найти передаточную функцию интегратора, полагая x(t) = et,
e – t, sin t, t.
7. Выполняя интегрирование, найти изображения по Лапласу
сигналов et, e – t, e – t + e t, 1 – e – t.
8. Найти реакцию апериодического звена на сигналы из п. 7 и
нарисовать соответствующие графики.
13
5. Варианты заданий
14
№
1
2
3
4
5
6
7
8
9
10
k
0,2
0,3
0,4
0,5
0,21
0,31
0,41
0,51
0,22
0,32
а
0,1
0,2
0,3
0,4
0,5
0,6
0,7
0,8
0,9
1,0
b
3,0
2,5
2,0
1,5
1,0
0,5
1,0
1,5
2,0
2,5
№
11
12
13
14
15
16
17
18
19
20
k
0,42
0,52
0,23
0,33
0,43
0,53
0,24
0,34
0,44
0,54
а
1,2
1,4
1,6
1,8
2,0
2,2
2,4
2,6
2,8
3,0
b
3,0
3,5
4,0
4,5
5,0
4,5
4,0
3,5
3,0
2,5
№
21
22
23
24
25
26
27
28
29
30
k
0,25
0,35
0,45
0,55
0,26
0,36
0,46
0,56
0,27
0,37
а
3,2
3,4
3,6
3,8
4,0
4,5
5,0
5,5
6,0
6,5
b
2,0
1,5
1,0
0,5
1,0
1,5
2,0
2,5
3,0
3,5
Лабораторная работа № 2
МОДЕЛИРОВАНИЕ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ
Цель работы: освоение методики моделирования линейных дифференциальных уравнений в системе MATLAB и SIMULINK.
1. Теоретические сведения
1.1. Линейное дифференциальное уравнение
Многие физические процессы, такие как колебания маятника,
движение стрелки гальванометра, изменение высоты при посадке
самолета, процессы в электрическом колебательном контуре могут
быть описаны линейным однородным дифференциальным уравнением второго порядка
(t) + a1x (t) + a0 x(t) =
x
0. (1)
Здесь а0, а1 – постоянные коэффициенты, определяющие характер
процесса, точкой обозначается производная по времени. Амплитуда
переменной x(t) зависит от начальных условий, например, от начального отклонения x0 маятника и его начальной скорости x 0 .
Вид теоретического решения дифференциального уравнения (1)
определяется корнями его характеристического полинома [1]
p2 + a1 p + a0 =
0.
Если корни вещественные и различные р1 = a1, р2 = a2, то решение имеет вид
=
x(t) C1 ea1t + C2 ea2t .
Если корни комплексные р1,2 = a ± ib, то решение имеет вид
=
x(t) C1 eat sin bt + C2 eat cos bt.
Постоянные С1 и С2 находят, подставляя начальные условия
в выражения для x(t) и x (t) при t = 0.
Пример 1. Дано дифференциальное уравнение
 + 2x + 2x = 0,
x
x0 = x 0 = 1.
Его характеристическое уравнение p2 + 2p + 2 = 0 имеет корни
p1,2 =-1 ± i. Следовательно, общее решение будет следующим:
=
x(t) C1 e-t sin t + C2 e-t cos t.
15
Дифференцируя, находим выражение для x :
x (t) =
-(C1 + C2 ) e-t sin t + (C1 - C2 ) e-t cos t.
При t = 0 с учетом начальных условий получаем C1 = 2, С2 = 1.
Следовательно,
=
x(t) e-t (2 sin t + cos t).
Эффективным средством решения дифференциальных уравнений является численное моделирование в одном из математических пакетов (MATHCAD, MATLAB, SIMULINK и др.). График
решения x(t) наблюдается на экране дисплея. В пакете MATLAB
для этой цели имеются команды initial, lsim, ode23, ode45, dsolve.
Дополнительныe возможности для пользователя предоставляет моделирование в SIMULINK.
1.2. Структурное моделирование
линейных дифференциальных уравнений
При структурном моделировании дифференциальных уравнений в пакете SIMULINK необходимо составить схему моделирования. На ней изображаются вычислительные блоки (усилители,
сумматоры, интеграторы) и связи между ними. При проведении моделирования эта схема набирается на экране дисплея с помощью
мыши или клавиатуры. По своему смыслу этот процесс аналогичен
вводу программы, однако он более прост и нагляден. Подробная информация о реализации таких схем в SIMULINK имеется в пособии
[2, разд. 3].
Рассмотрим методику составления схемы моделирования на
примере однородного линейного дифференциального уравнения
второго порядка
 + 2x + 3=
x
x 0, x(=
0) 2, x (=
0) 4. (2)
Для построения схемы моделирования воспользуемся методом
понижения производной (методом Кельвина). В нем можно выделить четыре шага.
Шаг 1. Разрешаем исходное уравнение относительно старшей
 =
производной. В частности для уравнения (2) получаем x
-2x - 3x .
Шаг 2. Полагаем старшую производную известной и выполняем
ее последовательное интегрирование, получая все низшие производные и саму переменную х. В случае уравнения (2) для этого потре16
4

x
2
x
x
–2
–3
Рис. 1. Схема моделирования уравнения (2)
буется два последовательно включенных интегратора, на выходах
которых получим сигналы x и x.
Шаг 3. Формируем старшую производную, используя уравнение,
полученное на первом шаге. В нашем примере для этого потребуется сумматор, складывающий сигналы x и x, домноженные, соответственно, на коэффициенты –2 и –3.
Шаг 4. Объединяем схемы, полученные на втором и третьем шагах, в общую схему моделирования, указываем начальные условия
интеграторов.
Применение этой методики для уравнения (2) приводит к схеме,
показанной на рис. 1. Она содержит два интегратора, два масштабных усилителя и сумматор (обозначен кружочком).
Выходной сигнал схемы подается на имитатор осциллографа
(блок Scope) или передается в рабочее пространство MATLAB (блоки OUT или ToWorkspaсe).
1.3. Системы линейных дифференциальных
уравнений первого порядка
Многие технические объекты можно описать системой n линейных дифференциальных уравнений первого порядка:
dY
= AY + bu, dt
(3)
где и – входной сигнал; Y – вектор-столбец выходных переменных
yi; b – вектор-столбец коэффициентов bi; A – квадратная матрица
коэффициентов aij, i, j = 1,n.
17
Например, при моделировании летательного аппарата составляющими вектора Y могут быть текущие координаты самолета
и скорости их изменения, тогда матрица A будет характеризовать
динамику самолета, а слагаемое bи описывать управляющие воздействия, формируемые летчиком или автопилотом.
Один из методов решения системы дифференциальных уравнений основан на предварительном переходе от системы (3) к одному
уравнению n-го порядка [1]. Для этого из уравнений системы и из
уравнений, полученных их дифференцированием, исключают все
переменные кроме одной. Для нее получают одно дифференциальное уравнение. Решая его, определяют эту переменную, а остальные
находят, по возможности, без интегрирования.
Пример 2. Дана система из двух дифференциальных уравнений
y1 = 3y1 - 2y2 + 5,
y2 =
-6y1 - y2 .
(4)
После дифференцирования первого уравнения получаем:
y1 = 3y1 - 2y2 = 21y1 - 4y2 + 15.
Чтобы исключить у2, вычтем отсюда удвоенное первое уравнение
системы (4):
y1 - 2y1 - 15y1 =
5.
Мы получили линейное неоднородное дифференциальное уравнение второго порядка. Общее решение этого уравнения представляет собой сумму общего решения соответствующего однородного
уравнения и частного решения=
y1 y1îäí + y1÷àñò . Так как корни характеристического уравнения р2 – 2р–15 = 0 вещественны и различны: р1 = –3, р2 = 5, то решение y1îäí имеет вид
1e
-3t
+ 2 e5t . Скла-
5
1
y1÷àñò =
=
- , получаем
15
3
y1 = C1e-3t + C2 e5t - 1 / 3. Переменную y2 находим из соотношения
=
y2 0,5(3y1 + 5 =
y1 ) 3Ñ1e-3t - Ñ2 e5t + 2.
Для определения постоянных коэффициентов С1 и С2 используют начальные условия системы. Аналогичным образом этот метод
применяется и для систем уравнений более высоких порядков.
дывая его с частным решением
18
1.4. Моделирование системы линейных
дифференциальных уравнений
Если задача описывается системой дифференциальных уравнений пер­вого порядка, то для ее моделирования по методу понижения
производной достаточно составить схемы для каждого уравнения отдельно и соединить их между собой. Например, схема моделирования системы уравнений (4) будет иметь вид, показанный на рис. 2.
Для наблюдения графиков сигналов у1(t), у2(t) в SIMULINK используется блок осциллографа SCOPE, а для наблюдения фазовой
траектории у2 = f (у1) – блок осциллографа XY Graph.
2. Задание по работе и содержание отчета
1. Теоретическое решение уравнения (1) при заданных значениях а0, а1 и начальных условиях x(0) = 5, x (0) = 0. Таблицы расчетных данных, графики решений x(t), x (t), график фазового портрета
x = f (x).
2. Схема моделирования заданного уравнения применительно
к SIMULINK.
3. Теоретическое решение системы дифференциальных уравнений (3) для случая
a11 a12 
1 
 y1 
0
, b =
, Y  ,=
, x 5, (5)
=
A 
=
Y (0)  =


0 a22 
0
20 
y2 




при заданных значениях аij. Графики решений у1(t) и у2(t) и график
фазового портрета у2 = f(y1).
Схема моделирования исходной системы уравнений применительно к SIMULINK.
y1
5
+
–6
y2
–
3
–2
Рис. 2. Схема моделирования системы уравнений (4)
19
3. Порядок выполнения лабораторной работы
1. Набрать в SIMULINK схему моделирования уравнения (1),
установить коэффициенты и начальные условия.
2. Получить осциллограммы x(t), x (t) и x = f (x), сравнить их
с теоретическими графиками. Варьировать шаг и метод интегрирования.
3. Набрать схему моделирования системы уравнений (3), установить коэффициенты и начальные условия (5).
4. Получить осциллограммы у1(t), у2(t) и у2 = f(y1), сравнить их
с теоретическими графиками. Варьировать шаг и метод интегрирования.
5. Выполнить моделирование системы уравнений (3) в MATLAB,
используя команду lsim. Cравнить графики, полученные в MATLAB
и SIMULINK.
4. Контрольные вопросы
1. Решить следующие линейные дифференциальные уравнения:
 + x =
(0) =
а) 
x - 2x
4,
x(0) =
1, x (0) =
2,
x
-2.
 + 4x = sin 2t,
б) x
x(0) = x (0) = 0.
 + 4x = cos 2t,
в) x
x(0) = x (0) = 0.
2. При каком значении а и при каких начальных условиях реше + àx + x =
ние уравнения x
0 имеет вид:
а) x = sin t; б) x = cos t; в) x = et;
г) x = e–t; д) x = e3t; е) x = et/3.
3. В чем заключается метод понижения производной? Пользуясь
этим методом, составить схемы моделирования для всех вариантов п. 1.
4. Используя метод понижения производной, составить схемы
моделирования следующих дифференциальных уравнений:
 -=
 = x , x(0) = x (0) = -1;
а) x
1 0, x(0=
) x (0=
) 1,5; б) x
 + 7x =
(0) =
в) 
x + 2x
-8, x(0) =
-10, x (0) =
-5, x
3;
(6)
=
г) x=
0, x=
(0) x=
(0) x
(0) 
x=
(0) 0, x(4)=
(0) -2.
5. Схема моделирования представляет собой кольцо из трех интеграторов с единичными коэффициентами и одинаковыми начальными условиями. Найти моделируемое дифференциальное уравнение и его аналитическое решение.
6. Как изменятся графики решения линейного однородного дифференциального уравнения при замене знаков всех начальных условий на противоположные?
20
7. Описать процедуру перехода от системы дифференциальных
уравнений к одному уравнению и обратную процедуру, рассмотрев
случай n = 3. Привести пример.
8. Составить схему моделирования и найти решение системы линейных дифференциальных уравнений
=
Y AY
=
, yi (0) 1, если матрица A имеет вид
0 1 1 
 0 1 1
1 1




а) A = 1 0 1  , б) A =  -1 0 1  , в) A = 
.
1 1

1 1 0 
 -1 -1 0 
5. Варианты заданий
№
1
2
3
4
5
6
7
8
9
10
11
12
a1
0,1
0,1
0,5
0,1
0,1
0,5
0,1
0,1
0,5
0,1
0,1
0,6
a0
0,4
1,6
4,8
0,5
1,8
5,0
0,6
2,0
5,4
0,7
2,2
5,8
a11 –1,0 –1,0 –1,0 –1,0 –1,0 –1,0 –1,0 –0,9 –0,9 –0,9 –0,9 –0,9
a12
1,0
0,8
0,7
0,6 0,57 0,4 0.35 1,0
0,8
0,7
0,6
0,5
a22 –2,0 –1,8 –1,7 –1,6 –1,5 –1,4 –1,3 –1,9 –1,7 –1,6 –1,5 –1,4
№
13
14
15
16
17
18
19
20
21
22
23
24
a1
a0
0,1
0,8
0,3
2,4
0
6,0
0,9
8,8
0,1
0,9
0,3
2,6
0,7
6,4
1.1
9,0
0,2
1,0
0,3
2,8
0,8
6,8
0,6
5,8
a11 –0,9 –0,9 –0,8 –0,8 –0,8 –0,8 –0,8 –0,8 –0,8 –0,5 –0,5 –0,5
a12
0,4
0,3
1,0
0,8
0,7
0,6
0.5
0,4
0,3
1,0
0,8
0,7
a22 –1,3 –1,6 –1,6 –1,6 –1,5 –1,4 –1,3 –1,2 –1,1 –1,5 –1,3 –1,2
21
Лабораторная работа №3
МОДЕЛИРОВАНИЕ СЛЕДЯЩЕЙ СИСТЕМЫ
Цель работы: освоить различные способы описания линейных динамических систем и методы их моделирования в пакетах
MATLAB и SIMULINK.
1. Теоретические сведения
1.1. Описание систем структурной схемой
и передаточной функцией
В инженерной практике используются различные способы задания динамических систем – с помощью структурных схем, передаточных функций, дифференциальных уравнений, а также с помощью частотных и временных характеристик. Проиллюстрируем их на примере следящей системы, структурная схема которой
представлена на рис. 1. В ее состав входят инерционное усилительное звено с передаточной функцией k/(T1p + 1), двигатель с передаточной функцией 1/(T2p) и вычитающее устройство для сравнения
входного сигнала u и выходного сигнала следящей системы y. Следящая система должна работать таким образом, чтобы угол поворота двигателя у по возможности точно равнялся значению входного
сигнала и (задача слежения).
Способ задания моделей объектов с помощью схемы (типа приведенной на рис.1) называется структурным, поскольку он отражает
реальную структуру объекта.
По передаточным функциям отдельных блоков можно построить
Y ( p)
, свяобщую передаточную функцию следящей системы Q( p) =
U ( p)
зывающую изображения по Лапласу входного и выходного сигналов.
u
+
e
k
T1 p + 1
x1
1
T2 p
–
Рис. 1. Структура следящей системы
22
x2
y
Для этого в соответствии со структурной схемой выписывается
система уравнений
1
k
Y ( p) = ⋅
⋅ e( p), e=
( p) U ( p) - Y ( p), (1)
T2 p (T1 p + 1)
которая затем преобразуется к одному уравнению, путем исключения переменной e(p):
k
=
Y ( p)
⋅ (U ( p) - Y ( p)).
T2 p ⋅ (T1 p + 1)
Далее, выражая выходной сигнал Y(p) через входной U(p), получаем
k
Y (p) =
⋅ U (p) = Q(p) ⋅ U (p).
T2 p ⋅ (T1 p + 1) + k
где Q(p) - передаточная функция системы.
В нашем случае она имеет вид
k
Q( p) =
.
2
T1T2 p + T2 p + k
По сравнению со структурным описанием передаточная функция является более компактной математической моделью, в то же
время она позволяет анализировать такие характеристики, как
устойчивость и минимальность объекта.
1.2. Описание систем дифференциальными уравнениями
От передаточной функции легко осуществить переход к описанию системы с помощью дифференциального уравнения. В рассматриваемом случае для этого достаточно в уравнении
(T1T2 p2 + T2 p + k)Y ( p) =
kU ( p)
раскрыть скобки и заменить оператор p оператором дифференцирования d/dt
T1T2 y + T2 y + ky =
ku. (2)
Решая это дифференциальное уравнение, можно найти реакцию
следящей системы на любое входное воздействие.
Аналитическое решение y(t) дифференциального уравнения (2)
является суммой решения однородного уравнения yодн(t) и частного
решения дифференциального уравнения yчастн(t).
23
Для получения yодн(t) составляем характеристическое уравнение
T1T2p2 + T2p + k = 0 и находим его корни p1 и p2. Если они вещественные и различные, то решение однородного уравнения ищется
в виде yîäí
=
(t) C1e p1t + C2e p2t , где С1 и С2 - коэффициенты, зависящие от начальных условий и определяемые в дальнейшем. Если
корни одинаковые (кратные) p1 = p2, то решение имеет вид
yîäí (=
t) (C1 + C2t)e p1t . Паре комплексных корней p1,2 = a ± jb соответствует решение =
yîäí (t) eat (C1 sin bt + C2 cos bt).
Во всех случаях система оказывается устойчивой, если корни лежат в левой полуплоскости (при этом решение однородного уравнения с течением времени стремится к нулю).
Частное решение дифференциального уравнения определяется
видом правой части дифференциального уравнения (2). Если, например, там стоит экспоненциальная функция u = e–t, то и частное
решение нужно искать в виде экспоненты yчастн = Ce–t. Если u = 1(t),
его следует искать в виде константы yчастн = C. Для определения C
надо подставить частное решение в дифференциальное уравнение.
Учитывая, что производная от константы равна нулю, находим, что
в последнем случае C = 1.
Значения постоянных С1, С2 определяются путем подстановки
в полученное решение начальных условий. Например, в случае нулевых начальных условий и решения вида y(t) = C1 e p1t + C2 e p2t + 1
постоянные С1 и С2 находятся из системы уравнений
C1 + C2 + 1 = 0; p1C1 + p2C2 = 0.
Наряду с заданием объекта одним дифференциальным уравнением типа (2) часто используют описание с помощью системы дифференциальных уравнений первого порядка. Оно известно как матричное описание или описание в пространстве состояний.
Для получения описания следящей системы в пространстве состояний выберем в качестве переменных состояния x1 и x2 выходные сигналы звеньев первого порядка на структурной схеме рис. 1.
Составим для каждого из них дифференциальное уравнение первого порядка
1
k
k
1
x1 =
- x1 - x2 + u, x2 =
x1..
T1
T1
T1
T2
Кроме того, запишем алгебраическое уравнение для выходного
сигнала y = x2.
24
В матричном виде это описание имеет вид
 AX + bu, y = cX,
=
X
где
 x1 
 -1 / T1
X =  , A = 
 1 / T2
 x2 
k / T1 
-k / T1 
, b=
 , c = 0 1 .

0 
 0 
Анализируя это описание, можно оценить устойчивость, управляемость, наблюдаемость и другие характеристики системы.
1.3. Взаимосвязь описаний
Все рассматриваемые виды описаний тесно взаимосвязаны, поэтому, зная одно из них, можно получить остальные. Например,
связь между матрицами A, b, c описания в пространстве состояний
и передаточной функцией системы Q(p) задается уравнением
Q=
( p) c( pE - A)-1 b, (3)
где p - оператор Лапласа, E - единичная матрица.
Любое из рассмотренных описаний системы позволяет рассчитывать ее реакцию на типовые входные сигналы. Чаще всего систему
характеризуют реакцией на дельта-функцию u = d(t) и на единичную функцию (функцию единичного скачка) u = 1(t). Эти реакции
известны как импульсная весовая характеристика системы y = q(t)
и переходная характеристика y = p(t). Их изображения по Лапласу
связаны с передаточной функцией формулами
q(t) ⇔ Q( p);
p(t) ⇔ Q( p)/ p, (4)
которые удобно использовать для нахождения весовой и переходной
характеристик.
Другой подход к описанию системы связан с использованием
частотных характеристик. Они получаются рассмотрением функции комплексной переменной, получаемой из формулы (3) заменой
p = jw: Q(jw) = c(jwE – A)–1b.
1.4. Моделирование в пакете MATLAB и SIMULINK
Пакет MATLAB поддерживает все виды описаний динамических систем, включая структурные схемы, передаточные функции и матричное описание в пространстве состояний. Для работы
со структурными схемами в пакете MATLAB имеется приложе25
ние SIMULINK. Его можно вызвать, набирая в командном окне
MATLAB команду simulink.
Численное моделирование следящей системы в MATLAB выполняется с помощью команд impulse, step, lsim. Предварительно надо
ввести числитель num и знаменатель den передаточной функции
либо матрицы A, B, C, D описания в пространстве состояний и сформировать структуру sys=tf(num,den) либо sys=ss(A,B,C,D). После
этого весовая функция и переходная функция находятся командами impulse(sys), step(sys), а реакции на произвольные входные сигналы, такие как u = e–t, рассчитываются с помощью команды lsim.
Реализация различных соединений блоков может быть осуществлена программно с помощью команд parallel, series, feedback,
append и некоторых других. Для этой цели можно использовать
также команды + , –, *.
В MATLAB можно получать не только численное, но и символьное
решение дифференциальных уравнений. Это делается с помощью команды dsolve тулбокса SYMBOLIC. Входными аргументами команды
служат дифференциальное уравнение и начальные условия. Напри2
мер, для решения дифференциального уравнения y + 2y + 3y =
с нулевыми начальными условиями следует набрать код
>> y=dsolve(‘D2у+3*Dу+2*y=2’, ‘Dy(0)=0’,’y(0)=0’)
MATLAB выдаст ответ y=1+exp(-2*t)-2*exp(-t), т. е. y = 1 + e–2t – 2e–t.
2. Задание по работе и содержание отчета
В работе исследуется динамика следящей системы, заданной
структурной схемой (рис. 1) при значениях параметров k, T1, T2, приведенных в таблице вариантов заданий.
Отчет по работе должен содержать:
1. Исходную схему моделирования с заданными численными
значениями параметров и передаточную функцию Q(p), полученную из уравнения (1).
2. Дифференциальное уравнение второго порядка, описывающее следящую систему, его аналитическое решение и график выходного сигнала y(t) при входном сигнале u = e–t и нулевых начальных условиях.
26
3. Описание следящей системы в пространстве состояний, передаточную функцию системы, полученную по формуле (3). Формулы
и графики весовой и переходной характеристик.
4. Схемы моделирования следящей системы применительно
к SIMULINK, содержащие осциллографы и генераторы входных
сигналов (для генерирования сигнала u = e–t использовать интегратор с обратной связью). Программы численного и символьного моделирования в MATLAB.
3. Порядок выполнения работы
1. С помощью пакета SIMULINK построить схему моделирования в соответствии с рис. 1.
Получить графики выходных сигналов (весовой функции, переходной функции и реакции на u = e–t). Проверить их соответствие
расчетным.
2. Параллельно со схемой моделирования следящей системы набрать модель передаточной функции Q(p) следящей системы и сравнить их выходные сигналы.
3. Выполнить моделирование в пакете MATLAB, используя разные
описания системы. Сравнить результаты моделирования в MATLAB
и SIMULINK.
4. Построить графики фазовых траекторий в плоскости (x1, x2)
(в SIMULINK для этого потребуется блок XY-graph). Построить частотные характеристики следящей системы, используя команды
bode, nyquist, ltiview.
4. Контрольные вопросы
1. Описать общую процедуру перехода от произвольной структурной схемы к системе линейных дифференциальных уравнений
первого порядка.
2. Найти реакцию своего варианта следящей системы на входной
сигнал u = t и построить график выходного сигнала.
3. Найти передаточную функцию следящей системы, если
устройство сравнения реализовано в соответствии с одной из следующих формул:
e= u - y, e= u - ky, e + e = u + u - y.
27
4. Как повлияет изменение знака обратной связи в следящей системе на ее устойчивость и вид переходной характеристики?
5. Найти передаточную функцию следящей системы, если передаточная функция двигателя равна
1
1
p
а)
,
б)
,
в)
.
T1 p + 1
T1 p - 1
T2 p + 1
6. Найти матрицы описания в пространстве состояния для пп. 3
и 5.
7. Сравнить графики весовой и переходной функций разомкнутой и замкнутой системы для своего варианта заданий.
5. Варианты заданий
№
1
2
3
4
5
6
7
3
15
k
4
16
21
T1
0,2
0,1
0,1
T2
5
20
5
4
8
№
11
12
13
14
k
6
8
10
T1
0,75
1,25
T2
4
4
28
7
8
9
10
12
5
10
16
0,2
0,4
0,8
1,2
4
8
2
2.5
5
15
16
17
18
19
20
12
14
16
20
18
16
12
1,5
2
2,5
2
1,5
0,8
0,4
0,8
2
3
2
5
10
5
8
5
0,125 0,25 0,125
Лабораторная работа № 4
МОДЕЛИРОВАНИЕ ДВОЙНОГО МАЯТНИКА
Цель работы: исследование колебаний двойного маятника путем
компьютерного моделирования в пакетах SIMULINK и MATLAB.
1. Теоретические сведения
1.1. Вывод дифференциальных уравнений
двойного маятника.
Динамика механических систем описывается законами Ньютона. При отсутствии трения они приводят к системе линейных дифференциальных уравнений второго порядка, матричная запись которых имеет вид
 + KX = 0, (1)
MX
где М, K – постоянные матрицы размеров m×m.
Если, например, речь идет о движении системы материальных тел,
 характеризуют положение и ускорение этих тел,
то векторы Х и X
а матрицы М и K зависят от масс и сил. Общее решение системы (1)
содержит n = 2m произвольных постоянных, для их определения
необходимо знать n начальных условий. Вид решения определяется
корнями p1, ..., pn характеристического уравнения, которое полу0.
чают, приравнивая нулю определитель Mp2 + K =
Для консервативных колебательных систем (систем без трения и
потерь энергии) корни получаются чисто мнимыми pi = ±jω1. Положительные числа w1, ..., wm называются циклическими частотами
системы.
Таким корням соответствует общее решение вида
X(t) = с1H1cos(w1t + j1) + ... + сmHmcos(wmt + jm),
(2)
где сi, ji – произвольные постоянные. Векторы Hi удовлетворяют алгебраическим уравнениям вида (Mpi2 + K)Hi =
0, которые получаются в результате подстановки отдельных компонент решения (2)
в систему (1).
В качестве примера механической системы, описываемой уравнениями вида (1), рассмотрим двойной маятник. Так называют систему из двух маятников – тяжелого, массы m1 и легкого, массы m2.
Легкий маятник подвешен к тяжелому, как это показано на рис. 1.
29
Массы m1 и m2 будем считать точечными,
длины нитей – одинаковыми l1 = l2 = l,
углы a1 и a2 – малыми.
α1
l1
При указанных условиях для двойного
маятника характерно так называемое явx1
ление биений, сопровождаемое цикличеm1
ским обменом энергией между маятникаh1
ми. Внешне картина колебаний выглядит
довольно неожиданно: без видимых приα2
l2
чин один из маятников время от времени
самопроизвольно останавливается, а друx2
гой начинает интенсивно раскачиваться.
m2
Подобные колебания могут возникать при
h2
спуске на парашюте, при подъеме по веревочной лестнице и в других ситуациях.
Рис 1. Двойной маятник
Для вывода дифференциальных уравнений малых колебаний двойного маятника воспользуемся законом сохранения полной энергии Е, согласно которому
Е = Ек + Eп = const,
(3)
где Ек и Еп – кинетическая и потенциальная энергии.
Выражение для кинетической энергии малых колебаний имеет вид
m1x12 m2 x22 (4)
+
,
2
2
где х1 и х2 – горизонтальные отклонения тяжелого и легкого маятников от положения равновесия. Потенциальная энергия определяется вертикальным отклонением маятников h1 и h2:
Eê ≅
Eï =m1 gh1 + m2 gh2 ≅
g
m x2 + m2 x12 + m2 (x2 - x1 )2  . 2l  1 1
(5)
Подставляя выражения (4) и (5) в (3) и переходя к матричной форме записи, получаем
где
E=
1 Ò  1 Ò
Õ MÕ + Õ KÕ = const, 2
2
m1 0 
 x1 
g m1 + 2m2 -m2 
=
M =
, K
=
, X  ,



m2 
l  -m2
 0 m2 
 x2 
знак Т означает транспонирование.
30
(6)
Геометрически уравнение (6) задает некоторый эллипсоид в четырехмерном пространстве состояний с координатами x1, x2 , x1, x2 .
Размеры эллипсоида пропорциональны полной энергии маятника.
Каждому состоянию маятника соответствует определенная точка
эллипсоида, при колебаниях она перемещается по некоторой траектории, лежащей на его поверхности.
Чтобы найти уравнение такой траектории, продифференцируем
равенство (6) по времени:
T
 T MX
 + X
=
 T (MX
 +
X
KX X
=
KX) 0. (7)
При дифференцировании мы воспользовались следующей формулой:
d T
 T MY + Y T MY
 = 2Y
 T MY
.
Y MY = Y
dt
Равенство (7) должно выполняться при любых значениях первого сомножителя, поэтому второй сомножитель должен равняться
 + KX = 0.
нулю MX
Найдена система дифференциальных уравнений, описывающих
движение двойного маятника. Ее более подробная запись имеет вид
1  g m1 + 2m2
m1 0   x
 0 m  x
+ 
2   2  l  -m2

-m2   x1  0 
.
=
m2   x2  0 
Переходя к скалярным уравнениям и вводя обозначения m2 = m2/m1,
окончательно получаем
k2 = g/l,
1 + k2 (1 + 2m2 )x1 - m2k2 x2 = 0,
x
2 - k2 x1 + k2 x2 =
x
0.
(8)
1.2. Решение дифференциальных уравнений
двойного маятника
Уравнения (8) можно решить аналитически, так как они представляют собой линейные однородные дифференциальные уравнения с постоянными коэффициентами. Для упрощения дальнейших
выкладок положим k = 1 и будем считать m << 1, тогда система уравнений (8) принимает вид:
 1 -m2 
=
X
+ AX 0=
,
A 
.  -1 1 
(9)
31
Выпишем характеристическое уравнение этой системы
Ep2 + =
A
p2 + 1
-1
-m2
2
= ( p2 + 1)2 - m=
0.
p2 + 1
Оно имеет 4 чисто мнимых корня p1,2 = ±jw1, p3,4 = ±jw2, где циклические частоты w1, w2 определяются равенством ω1,=
1 ± m ≈ 1 ± 0,5m.
2
Следовательно, общее решение имеет вид (2) при m = 2:
X(t) = C1H1cos(w1t + j1) + C2H2cos(w2t + j2).
Векторы H1 и Н2 – это собственные векторы матрицы А, они имеют вид:
 1 
 1 
=
H1 =
,
H2 

.
 -1 / m 
1 / m 
Для определения произвольных постоянных с1, с2, j1, j2 зада=
дим начальные условия х1(0) = 1; х2(0) = 0, x=
1 (0) x
2 (0) 0.
 , получаем j = j = 0;
Подставляя их в выражение для Х и Õ
1 2 С1 = С2 = 0,5. Окончательный вид решения:
x1 = 0,5 (cosw1t + cosw2t) = cos(µt) cos t,
x2 = (– 0,5/µ)(cosw1t – cosw2t) = (1/µ)sin(µt) sin t.
(10)
В отчете требуется построить графики по формулам типа (10).
При малых m первый сомножитель в этих формулах меняется значительно медленнее, чем второй, что позволяет рассматривать его как
огибающую результирующего графика. Поэтому при построении
графика х1(t) удобно сначала построить огибающие ±cosµt, а затем
заполнить область между ними косинусоидальным сигналом с периодом 2p. Аналогично строится график функции х2(t).
Графики х1(t) и х2(t) имеют экстремумы и нули в точках, кратных p/2, поэтому целесообразно рассчитать их значения для этих
моментов и полученные данные свести в таблицу. Этой таблицей
удобно пользоваться и при построении фазовой траектории. Пример
графиков функций х1(t) и х2(t) для m = 0,125 показан на рис. 2.
Графики представляют собой «быстрые» колебания с периодом 2p,
модулированные «медленными колебаниями» с периодом 2pm, и наглядно описывают явление биений, заключающееся в циклической
«перекачке» энергии от одного маятника к другому.
Биения можно охарактеризовать тремя параметрами – периодом
быстрых колебаний t, периодом медленных колебаний Т и числом
«быстрых» колебаний за период биений n = T/t.
32
x1
x2
1
0,8
0,6
0,4
0,2
0
–0,2
–0,4
–0,6
–0,8
–1
10
8
6
4
2
0
–2
–4
–6
–8
–10
t
T
t
τ
Рис. 2. Графики колебаний двойного маятника
Построение графика фазовой траектории также удобно начинать
с нахождения его «огибающей», т. е. геометрической фигуры, внутри которой он расположен. Для двойного маятника такой фигурой
является ромб с центром в начале координат.
Чтобы построить фазовую траекторию, следует сначала нарисовать этот ромб, а затем, пользуясь полученными ранее таблицей
и графиками рис. 2, последовательно наносить точки (х1, х2), соответствующие экстремумам и нулям функций х1 и х2, и соединять
их плавной кривой, вписанной в ромб. Если отказаться от упрощающего предположения m<<1 то ромб перестанет быть строго симметричным относительно координатных осей и примет вид, показанный на рис. 3.
1.3. Составление структурной схемы моделирования
Для построения такой схемы воспользуемся методом понижения
производных, применяя его к каждому из уравнений (8). В соответствии с этим методом предположим, что нам известны вторые про1, x
2 . Пропуская каждую из них через цепочку из двух
изводные x
последовательно включенных интеграторов, получим переменные
x1 и x2. После этого сформируем необходимые нам значения вторых
производных на основе равенств (8):
1 = -k2 (1 + 2m2 )x1 + m2k2 x2 ,
x
2 = -k2 - x1k2 x2 .
x
33
x2 10
1
–1
x1
–10
Рис. 3. Траектория в плоскости x1, x2
1
и
–
k2
х4
х2
µ2k2
х3
х1
–
v2k2
Рис. 4. Схема моделирования двойного маятника
В результате получаем схему, показанную на рис. 4.
Таким образом, схема моделирования двойного маятника состоит из двух чисто колебательных звеньев с близкими собственными
частотами k и vk (v2 = 1 + 2µ2), соединенных в «кольцо». Левая схема моделирует колебания легкого маятника, правая – тяжелого,
взаимное влияние маятников учитывается связями между схемами (коэффициенты k2 и m2k2). Начальные условия, показанные на
схеме, означают, что в первый момент тяжелый маятник отклоняют
от положения равновесия и без толчка отпускают, легкий маятник
при этом имеет нулевые значения скорости и координаты.
34
Используя схему рис. 4, легко получить матрицы описания двойного маятника в пространстве состояний:
 0

 0
A= 2 2
-n k

 k2

0
0
m2 k2
-k2
1 0
0 1 
,
0 0

0 0 
1 0 0 0 
C= 
.
0 1 0 0 
Введение вектора B = [0 0 0 1]T позволяет рассматривать движение двойного маятника при наличии управляющего воздействия,
приложенного к нижнему маятнику.
2. Задание по работе и содержание отчета
В лабораторной работе осуществляется компьютерное моделирование двойного маятника в пакетах SIMULINK и MATLAB. Основой для моделирования является схема рис. 4 и описание в пространстве состояний. Численные значения параметров двойного маятника приведены в таблице вариантов.
Отчет должен содержать:
1. Теоретическое решение системы уравнений (8) при заданных
.
значениях k, m2 для начальных условий х1(0) = 5; х2(0) = х1(0) = .
= х2(0) = 0.
2. Расчет численных значений параметров Т, t, n, таблицу и графики функций х1(t), х2(t), график фазовой траектории х2 = f(х1);
графики должны отражать полтора – два периода биений.
3. Матрицы A, B, C описания системы в пространстве состояний
в случае приложения управляющего воздействия к нижнему маятнику и передаточную функцию системы для этого случая.
4. Схему моделирования исходной системы в SIMULINK и программу моделирования на языке MATLAB.
3. Порядок выполнения работы
1. Выполнить в пакете SIMULINK раздельное моделирование
колебательных звеньев (маятников) (рис. 4), не соединяя их между собой. Принимая n = 1, установить нулевые начальные условия.
Подать на входы обеих звеньев единичный входной сигнал и убедиться, что колебания на выходах звеньев совпадают (полупериод
колебаний должен равняться t = p/k).
35
2. Отключить единичный сигнал, соединить звенья между собой
«в кольцо», согласно рис. 4, установить начальные условия х1(0) = 1.
Наблюдать графики сигналов х1(t), х2(t) и их разности. Сравнить
экспериментальные оценки величин Т, t, n с их теоретическими
значениями.
3. Выполнить моделирование двойного маятника в пакете
MATLAB, используя описание в пространстве состояний и приняв
n = 1. Наблюдать фазовую траекторию х2 = f(х1) (на экране должен
быть виден ромб).
4. Установить коэффициент n2 = 1 + 2m2, наблюдать изменения
фазового портрета («перекос» ромба).
5. Найти начальные условия, при которых оба маятника качаются синхронно (синфазно и противофазно). Исследовать поведение
двойного маятника в случае приложения управляющего воздействия к нижнему маятнику.
4. Контрольные вопросы
1. Механическая система, содержащая две точечные массы m1,
m2 и три пружины (см. рис. 5), совершает колебания.
Требуется, используя закон сохранения энергии, найти дифференциальное уравнение малых колебаний, определить циклические
частоты и составить схему моделирования. Принять m1 = m2 = 1,
k = 1, трением пренебречь. Для определения потенциальной энерkx2
гии сжатой пружины использовать формулу Eï
.
2
m1
k3
m2
k2
k1
Рис. 5. Двухмассовая система с закрепленными концами
m1
k
m2
Рис. 6. Двухмассовая система со свободными концами
36
2. Две массы, соединенные пружиной, лежат на полированном
столе (рис. 6). Их прижимают друг к другу и отпускают. Составить
дифференциальное уравнение и исследовать движение такой системы, полагая m1 = m2 = m, k = 1.
Трением пренебречь.
3. Какой вид примет теоретическое решение (10) в случае n2 ≠ 1?
4. Найти вид теоретического решения и график фазового портрета уравнений (9) при следующих значениях начальных условий:
1
m 
m
а) X(0) =   ,
б) X(0) =   ,
в) X(0) =   .
 -1
1
1 
5. Какие условия надо выполнить, чтобы обеспечить колебания
двойного маятника с первой циклической частотой? Со второй циклической частотой?
6. При соединении схем моделирования, показанных на рис. 4,
между ними по ошибке включили инвертор. Найти вид сигнала x1
в этом случае, если m = 0,1, k = 1, n = 1.
7. Вывести уравнения (8) с помощью второго закона Ньютона.
5. Варианты заданий
№
k,
c–1
1
2
3
4
5
6
7
8
9
10
11
2
3
5
6
8
2
4
6
4
2
9
100m2
0,5
1,5
2,0
1,0
0,5
2,0
1,5
2,0
3,5
2,5
0,5
№
12
13
14
15
16
17
18
19
20
21
22
k, c–1
7
5
3
4
2
7
2
3
4
6
9
100m2
1,0
1,5
0,5
2,0
1,5
0,5
4,0
2,0
3,0
0,5
1,0
№
23
24
25
26
27
28
29
30
31
32
33
3
2
5
4
6
4
2
5
7
3
5
2,5
1,0
3,0
2,5
1,5
0,5
3,0
1,0
1,5
3,0
0,5
k,
c–1
100m2
37
Лабораторная работа № 5
МОДЕЛИРОВАНИЕ КОЛЕБАНИЙ СТРУНЫ
Цель работы: ознакомиться с методикой компьютерного моделирования колебаний струны с помощью пакетов MATLAB и
SIMULINK.
1. Теоретические сведения
1.2. Собственные частоты и собственные колебания струны
Задача о поперечных колебаниях струны – одна из классических
задач математической физики. В ней рассматривается туго натянутая струна, закрепленная в конечных точках. Если вывести струну из положения равновесия, например, оттянуть ее и затем отпустить, то струна начнет колебаться. В процессе колебания величина
отклонения u будет зависеть от абсциссы точки струны х и времени t. Таким образом, чтобы знать положение любой точки струны
в произвольный момент времени, надо найти зависимость u от х и t,
т. е. найти функцию u(x,t). Если считать струну абсолютно гибкой,
упругой (подчиняющейся закону Гука) и рассматривать малые колебания, то функция u(x,t) удовлетворяет уравнению:
a2
∂ 2u
∂ 2u
T
=
, a2 = , 2
2
r
∂x
∂t
(1)
где Т – сила натяжения струны, r – ее линейная плотность (масса
единицы ее длины).
Уравнение (1) называется уравнением свободных колебаний
струны или одномерным волновым уравнением. К нему сводится
не только рассматриваемая задача, но и многие другие, например,
задача о поперечных колебаниях летательных аппаратов, задача
о флаттере (колебаниях крыла самолета).
Для выделения конкретного решения в задаче о колебании струны задаются граничные условия двух видов: начальные и краевые.
Начальные условия показывают, в каком состоянии находилась
струна в момент начала колебания (т. е. описывают форму струны
и скорость ее точек при t = 0). Начальное состояние струны задается
двумя функциями:
∂u
u t =0 = f (x),
= j (x). (2)
∂ t t =0
38
Краевые условия показывают, что происходит на концах струны
во время колебаний. Если концы струны закреплены (начало струны – в начале координат, а конец – в точке х = l), функция u(x, t) будет подчиняться условиям:
=
u
0=
,
u
0. (3)
=
x 0=
x l
Характер колебаний струны сильно зависит от начальных условий. Наиболее простой случай получается, когда в начальный
момент времени струне придают форму полуволны синусоиды
f(x) = sinπx / l и отпускают без начальной скорости. Тогда все точки
струны будут совершать гармонические колебания с одной и той же
частотой ω1, так что в любой момент времени форма струны будет
отличаться от исходной только амплитудой (рис. 1, а).
Математически такие колебания описываются формулой
u ( x, t ) =p
A1 ( t ) sin x / l,
A1 ( t ) =
cos ω1t. (4)
Это так называемое первое собственное колебание струны, ему
соответствует низкочастотный чистый тон. Именно он используется музыкантами при настройке гитары или рояля. Чтобы найти
первую собственную частоту ω1, подставим функцию (4) в уравнение (1). После дифференцирования и сокращения на общий множитель получаем:
a2 p2 / l2 = ω12 ,
ω1 = ap / l.
Таким образом, решение уравнения (1) в рассматриваемом случае описывается формулой
u=
1 ( x, t ) cos
ap
p
t ⋅ sin x. l
l
(5)
Второе собственное колебание получим, задав начальное условие в виде f(x) = sin2πx/l. Форма струны в различные моменты вре-
u1
а)
u2
б)
x
x
0
l
0
l
Рис. 1. Первое и второе собственные колебания струны
39
мени для этого случая показана на рис. 1, б. Она описывается уравнением
u ( x, t ) =p
A2 ( t ) sin 2 x / l,
A2 ( t ) =
cos ω2t.
Подставляя эти выражения в уравнение (1), находим, что вторая
собственная частота ω2 = 2aπ/l. Отсюда получаем формулу для второго собственного колебания
=
u2 ( x, t ) cos
2ap
2p
t ⋅ sin x. l
l
(6)
Аналогично выводятся формулы для третьего и остальных собственных колебаний:
=
u3 ( x, t ) cos
3ap
3p
t ⋅ sin x, .
l
l
Если струна колеблется с собственной частотой ωk = kaπ/l, то
у нее будут k + 1 неподвижных точек (узлов) и k точек пучности (точек, где отклонения достигают максимума). Колебания такого вида
называются стоячими волнами.
Когда струна колеблется, она издает звук, высота которого возрастает с частотой колебаний. Если струна совершает собственные колебания, то самый низкий тон будет, когда частота равна w1.
Остальные тона, соответствующие частотам wk, называются обертонами или гармониками. Если струна совершает свободные колебания, то функция u(x, t), представляется в виде суммы отдельных гармоник. Это позволяет записать общее решение уравнения (1)
в виде бесконечного ряда Фурье.
1.2. Дискретизация пространственной координаты
Для перехода от уравнения (1) к конечномерной модели осуществим дискретизацию задачи по переменной x. В результате получится система обыкновенных дифференциальных уравнений, число уравнений в ней зависит от шага дискретизации h. С этой целью
выделим на струне n равноотстоящих точек с координатами х1, х2,
..., хn (рис. 2).
Тем самым струна условно разбивается на n равных участков длины h = l/n и заменяется приближенной моделью в виде n – 1 масс
(точек, бусинок), закрепленных на невесомой нити. Движение
рассматриваемых точек в процессе колебаний струны будет описываться некоторыми функциями времени u1(t), u2(t), ..., un–1(t).
40
u
...
x1
x2
x3 ... xi–1
...
x
xi ... xn–2 xn–1 xn = l
Рис. 2. Конечномерная модель струны
Для вычисления второй производной по х от функции и воспользуемся приближенной формулой
∂2u(x,t)
∂x
2
≈
x = xi
ui +1 - 2ui + ui -1
h
2
,
l
h=
.
n
Записывая теперь уравнение (1) для точек x1, x2, ..., xn–1, получаем систему обыкновенных дифференциальных уравнений второго
порядка
d2ui
a2
=
b
(
u
2
u
+
u
),
i
=
1
,
n
1
,
b
=
(7)
i +1
i
i -1
dt2
h2
с граничными условиями u0 = 0, un = 0, ui(0) = f(xi). Она представляет собой конечномерную модель струны.
В частности, для n = 4, b = 1 эта модель будет содержать три уравнения
1 =u2 - 2u1 + u0 ,
u

u2 =u3 - 2u2 + u1, (8)
3 =u4 - 2u3 + u2 ,
u
с краевыми условиями и0 = и4 = 0.
Суммарный порядок уравнений равен 6, т. е. мы имеем дело с
шестимерной моделью.
Начальную скорость струны будем считать нулевой
u=
(
1 0 ) u=
2 ( 0 ) u=
3 ( 0 ) 0. Чтобы найти решение этой системы, перепишем ее в матричной форме
 = A U,
U
0
 -2 1 0 
 u1 


 
A 0 =2 1 , U =
1
u2  .  0 1 -2
u3 
(9)
41
Решение будем искать в виде U = Hcosωt, где вектор Н и частота ω подлежат определению. После подстановки этого выражения
в (9) получаем соотношение
A 0 H = -ω2H.
Это означает, что Н и –ω2 – собственные векторы и собственные
числа матрицы А0. Чтобы найти собственные числа, выписываем
характеристический полином матрицы А0:
det ( lE - A0 ) = l3 + 6l2 + 10l + 4 = (l + 2)(l2 + 4l + 2),
и находим его корни
l1 =-2 + 2 =-0,5858, l2 =-2, l3 =-2 - 2 =-3,4142.
Им соответствуют собственные векторы, удовлетворяющие алгебраическим уравнениям A0Hi = liHi:
1 
 1 
 -1
 


 
H1 =  2  , H2 =  0  , H3 =  - 2  ,
 


 1 
1 
 1 
и собственные частоты ωi= -l i :
ω1 = 0,7654, ω2 = 1,4142, ω3 = 1,8478.
Общее решение системы (9) имеет вид
=
U ( t ) c1H1 cos ω1t + c2H2 cos ω2t + c3H3 cos ω3t, (10)
постоянные с1, с2, с3 зависят от начальных условий и1(0), и2(0), и3(0).
Уравнения для их определения получаем из соотношения (10)
при t = 0:
U ( 0 ) = c1H1 + c2H2 + c3H3 .
В нашем случае получаем систему
c1 - c2 + c3 =
u1 ( 0 ),
2c1 - 2c3 =
u2 ( 0 ),
c1 + c2 + c3 =
u3 ( 0 ).
Она легко решается:
1
u1 ( 0 ) + 2u2 ( 0 ) + u3 ( 0 )  ,
4
1
1
c2 = u3 ( 0 ) - u1 ( 0 )  , c3 = u1 ( 0 ) - 2u2 ( 0 ) + u3 ( 0 )  .
2
4
c1 =
42
Этим завершается получение аналитического решения для дискретной модели струны (8).
Особенно простые решения получаются, когда начальная форма
струны симметрична относительно середины u1(0) = u3(0), при этом
свойство симметрии сохраняется и в процессе колебаний. Например, при u1(0) = u3(0) = 2, u2(0) = 0 получим
u1 = u3 = cos ω1t + cos ω3t,
u2 =
2 cos ω1t - 2 cos ω3t.
1.3. Моделирование в SIMULINK
Для структурного моделирования системы уравнений (8) воспользуемся методом Кельвина, применяя его к каждому из трех
уравнений по отдельности и затем соединяя полученные схемы
между собой. В результате получаем схему моделирования на сумматорах и интеграторах, показанную на рис. 3. В его нижней части условно изображена исследуемая струна, разбитая на 4 участка.
Схема содержит три фрагмента, соединенных прямыми и обратными связями. Выход первого фрагмента u1 характеризует колебания
точки струны х1, выход второго фрагмента u2 характеризует колебания точки х2 и т. д.
Выходы внутренних интеграторов каждого из трех фрагментов
показывают, как изменяются скорости тех же точек.
При реализации схемы в SIMULINK нужно на входах интеграторов И1, И3, И5 поставить сумматоры, а выходные сигналы интеграторов И2, И4, И6 наблюдать с помощью блока осциллографа Scope.
Перед началом моделирования на интеграторах И2, И4, И6 устанавливаются начальные условия, соответствующие начальной форме
струны, а на интеграторах И1, И3, И5 – условия, соответствующие
начальным скоростям точек х1, х2, х3 (в нашем случае они равны
нулю).
u0 = 0
И1
И2
u1
–2
0
И3
И4
u2
–2
x1
u3
И6
И5
u4= 0
–2
x2
x3
l
Рис. 3. Схема моделирования колебаний струны
43
1.4. Моделирование в MATLAB
Для численного моделирования уравнений (8) в MATLAB удобнее всего использовать команду initial. Предварительно надо ввести
матрицы A, B, C, D описания в пространстве состояний и сформировать структуру sys = ss(A,B,C,D). В нашем случае эти матрицы имеют вид
0 0 0
0 0 0

0 0 0
=
A 
 -2 1 0
 1 -2 1

 0 1 -2
1
0
0
0
0
0
0
1
0
0
0
0
0
0 
1
=
, B
0
0

0 
0 
0 
 
0 
=
, C
0 
0 
 
0 
1 0 0 0 0 0 


=
0 1 0 0 0 0
, D
0 0 1 0 0 0 
0 
 
0  .
0 
Они могут быть получены из системы (8) после введения переменных u4, u5, u6, равных скоростям точек х1, х2, х3, либо по схеме рис. 3 путем выписывания уравнений для каждого интегратора.
Далее вводится массив времени t и вектор начальных условий U0
(его первые три компоненты – заданные числа, следующие три – нули). При формировании матриц в MATLAB можно пользоваться командами zeros, eye, ones, например:
>>B=zeros (6,1),
>>С=[eye(3), zeros(3)].
Результат выполнения команды
>>U=initial(sys,U0,t);
можно наблюдать с помощью команд plot(t, U) – обычные графики,
и mesh (U) – графики поверхности u(x, t) в трехмерном пространстве.
В MATLAB можно получить не только численное, но и символьное решение системы дифференциальных уравнений (8). Это делается с помощью команды dsolve тулбокса SYMBOLIC. Входными
аргументами команды служат дифференциальные уравнения и начальные условия (те и другие заключаются в одиночные кавычки –
строковый формат).
В нашем случае набираем код
>> y=dsolve(‘D2u1=u2-2*u1’,’D2u2=u3-2*u2+u1’,’D2u3=u2-2*u3’,’Du1(0)=
0’,’Du2(0)=0’,’Du3(0)=0’)
>> simplify(y.u1),simplify(y.u2),simplify(y.u3)
44
Результатом выполнения первой команды будет символьная
структура. Второй командой извлекаем из нее переменные u1, u2,
u3 и упрощаем их. Полученные формулы с точностью до обозначений совпадают с решением, найденным ранее в п. 1.2.
Для визуализации символьных решений служат команды ezplot
и ezsurf. В качестве примера их использования на рис. 4 показаны
поверхности, отвечающие первому и второму собственным колебаниям струны, описываемым формулами (5), (6).
Они построены с помощью команд
>>ezsurf(‘sin(pi*x)*cos(pi/4*t)’,[0,20,0,1])
>>ezsurf(‘sin(2*pi*x)*cos(pi/2*t)’,[0,20,0,1])
Погрешность, вносимую дискретизацией пространственной координаты, иллюстрирует рис. 5.
На нем показаны графики второго собственного колебания, полученные путем решения уравнения в частных производных (1)
1
и системы дифференциальных уравнений (8) при
=
a =
, n 4. Сплош4
sin(πx) cos(π/4t)
1
20
0
0
–1
1
0,8
0,6
0,4
0,2
t
x 0
sin(2πx) cos(π/2t)
1
0
–1
20
0,5
x
0
10
t
20
Рис. 4. Пространственное изображение собственных колебаний струны
45
1
0,5
cos(πt/2)
0
cos(1,41t)
–0,5
–1
0
2
4
6
8
10 t
Рис. 5. Графики собственного колебания для непрерывной
и дискретной моделей
ной кривой отвечает точное значение второй собственной частоты ω2 = π/2 = 1,57, пунктирной – значение ω2= 2= 1,41, равное
корню из второго собственного числа матрицы А0 (9). Для построения графика использовались команды
>>ezplot(‘cos(sqrt(2)*t)’,[0,10]);grid,
>>hold;ezplot(‘cos(pi/2*t)’,[0,10])
Очевидно, что с увеличением числа участков, на которые разбивается струна при дискретизации, кривые будут сближаться.
2. Задание по работе и содержание отчета
1. Изобразить на графике начальную форму струны для своего
варианта. Провести дискретизацию уравнения (1), условно разбив
струну на четыре равных участка.
2. Решить полученную систему уравнений теоретически. Для заданной f(x) построить графики ui(t), для трех точек струны х = 0,25;
0,5; 0,75.
3. Составить полную схему моделирования для SIMULINK. Рассчитать амплитуду и период косинусоиды на выходе каждого из
трех фрагментов схемы при отсутствии связей между ними.
4. Рассчитать начальные условия U1(0) и U2(0) для получения
первого и второго собственных колебаний шестимерной модели
струны. Оценить, на сколько процентов отличаются их собственные
частоты от теоретических.
46
5. Выписать матрицы А, В, С описания в пространстве состояний, привести программы численного и символьного моделирования в MATLAB.
3. Порядок выполнения лабораторной работы
1. Собрать в SIMULINK три фрагмента схемы моделирования, не
соединяя их между собой. Проверить их идентичность, установив
везде одинаковые начальные условия и сравнить выходы схем.
2. Соединить схемы между собой, установить начальные условия
для получения первого и второго собственных колебаний. Сравнить
их графики с теоретическими.
3. Установить начальные условия для своего варианта, полученные графики сравнить с теоретическими.
4. Выполнить моделирование в MATLAB, используя матрицы
описания в пространстве состояний и команду initial. Построить
график поверхности u(x, t) в трехмерном пространстве, используя
команду mesh.
5. Получить символьное решение в MATLAB, используя команду dsolve. Визуализировать полученное решение с помощью команды ezsurf.
4. Контрольные вопросы
1. Чем определяется общее число граничных условий? Чем краевые условия отличаются от начальных? Какой вид имеют начальные условия в задаче о колебаниях струны при игре на рояле и на
гитаре?
2. Вывести формулы конечно-разностной аппроксимации первой и второй производных. Оценить точность этих формул на конкретном примере.
3. Что такое собственные колебания и собственные частоты? Какой физический смысл этих понятий? Нарисуйте для различных
моментов времени форму струны, колеблющейся с третьей собственной частотой.
4. Выполнить дискретизацию, разбив струну на 3 участка и приняв b = 1. Для полученной системы дифференциальных уравнений: а) найти собственные частоты колебаний; б) найти собственные векторы матрицы А0; в) найти начальные условия для получения собственных колебаний; г) найти решение для случая
=
u1 ( 0 ) 0=
, u2 ( 0 ) 1.
47
5. Найти решение и нарисовать графики колебаний струны для
своего варианта при u1 ( 0 ) = 1, u2 ( 0 ) = 0, u3 ( 0 ) = -1.
6. Провести для уравнения струны дискретизацию не по пространственной переменной х, а по времени t. Сравнить полученную
схему моделирования со схемой для дискретизации по пространству. Указать физический смысл начальных условий для четных и
нечетных интеграторов. Какие сигналы следует подавать на вход
первого и последнего фрагментов схемы?
5. Варианты заданий
Рассматривается струна единичной длины (l = 1) с закрепленными концами. Ее движение описывается волновым уравнением (1)
с нулевыми граничными условиями (3). Начальные условия задаются формулами (2), причем
j(x) =
0, f (=
x) a1 sin px + a2 sin 3p x, 0 ≤ x ≤ l.
Значения коэффициентов а1, а2 и параметра b, входящего в уравнения (7), приведены в таблице.
48
№
1
2
3
4
5
6
7
8
9
10
а1
1
1
0,5
1
0,5
1
1
0,2
1
–1
а2
1
0,5
1
1
0,5
1
1
–1
0,1
1
b
1
1
1
2
2
2
3
3
3
3
№
11
12
13
14
15
16
17
18
19
20
а1
0,7
0,9
–1
0,5
0,6
2
0,3
1
0,4
1
а2
0,8
1
–1
1
–0,5
0,1
1
0,3
1
0,8
b
1
1
0,5
0,5
0,5
0,5
1
1
1
2
№
21
22
23
24
25
26
27
28
29
30
а1
0,7
0,7
–0,7
0,5
–0,5
0,4
0,8
0,8
0,8
0,6
а2
0,8
–0,8
0,8
0,5
0,6
0,4
0,3
–0,1
–0,3
0,7
b
2
2
1
4
4
4
4
2
2
2
Лабораторная работа № 6
АНАЛИЗ СИСТЕМЫ УПРАВЛЕНИЯ
Цель работы: исследовать теоретически и с помощью пакета MATLAB устойчивость, управляемость, наблюдаемость и минимальность системы управления, заданной структурной схемой, и построить ее модель пониженного порядка.
1. Теоретические сведения [1]
1.1 Анализ устойчивости
Понятие устойчивости является одним из основных в теории
управления. Устойчивой называют систему, которая, будучи выведена из состояния равновесия, стремится вновь вернуться в это
состояние. Для неустойчивых систем характерна обратная тенденция. Примером неустойчивой системы может быть карандаш, стоящий на острие. Обычно исследование устойчивости систем сводится к анализу устойчивости соответствующих дифференциальных
уравнений. У неустойчивых уравнений решение неограниченно
возрастает со временем. В теории автоматического управления существуют различные методы анализа устойчивости. Для линейных
систем разработаны критерии устойчивости, которые можно разделить на корневые, алгебраические и частотные.
Корневой критерий устойчивости. Для того чтобы линейная динамическая система была устойчивой, необходимо и достаточно,
чтобы все корни ее характеристического уравнения лежали в левой
комплексной полуплоскости.
Например, математическая модель системы, описываемой дифференциальным уравнением
 + 2x + 5x =
x
0,
устойчива, так как корни характеристического уравнения
р2 + 2р + 5 = 0 имеют отрицательные вещественные части
p1,2 = –1±2i, т. е. лежат слева от оси ординат.
Алгебраический критерий устойчивости. Алгебраические критерии устойчивости позволяют судить об устойчивости на основе
анализа коэффициентов характеристического полинома. Наибольшую известность получил критерий устойчивости Гурвица. Из него, в частности, следует, что все коэффициенты устойчивого дифференциального уравнения должны быть положительны.
49
Для уравнений второго порядка это необходимое и достаточное
условие.
Для уравнений третьего порядка
 + a1x + a0 x =
a3
x + a2 x
0,
помимо положительности коэффициентов, должно выполняться дополнительное условие, а именно, произведение средних коэффициентов должно быть больше произведения крайних: а2а1 > а0а3.
 + 0,5x + 3x =
Например, уравнение 
x + 2x
0 неустойчиво, так как
2 × 0,5 < 1 × 3.
Частотный критерий устойчивости. Частотные критерии носят
графический характер. Они опираются на анализ графиков частотных характеристик – АЧХ, ФЧХ, АФХ (последняя известна также
как диаграмма Найквиста). Часть из них позволяет делать заключение об устойчивости замкнутой системы управления по частотным характеристикам разомкнутой системы.
Приведем в качестве примера критерий Найквиста. Обозначим
передаточную функцию разомкнутой системы Q(p) и охватим ее
единичной отрицательной обратной связью (рис. 1).
Пусть известно, что разомкнутая система устойчива. Тогда для
устойчивости замкнутой системы необходимо и достаточно, чтобы
годограф Найквиста разомкнутой системы не охватывал точку с координатами (–1, 0j) на комплексной плоскости.
В тулбоксе CONTROL пакета MATLAB диаграмма Найквиста системы sys строится командой nyquist(sys), указанная точка помечена
на ней красным крестиком.
В лабораторной работе исследуется система управления, заданная структурной схемой, поэтому для анализа устойчивости надо сначала найти ее передаточную функцию. Структурная схема исследуемой системы приведена на рис. 2. В ее состав входят
три апериодических звена, а также суммирующие и вычитающие
звенья.
В схеме можно указать три пути от входа до выхода с передаточными функциями
k1
-k2
k3
-k1
ak3
Q1 ( p) =
, Q2 ( p) =⋅
, Q3 ( p) =⋅
.
T1 p + 1
T2 p + 1 T3 p + 1
T1 p + 1 T3 p + 1
Общая передаточная функция получается, как их сумма
50
Q( p) = Q1 ( p) + Q2 ( p) + Q3 ( p) =
B( p)
.
A ( p)
(1)
u
y
Q(p)
Рис. 1. Система с обратной связью
k1
T1 p+ 1
x1
y
–
u
a
+
k2
T2 p+ 1
x2
+
k3
T3 p+ 1
x3
Рис. 2. Структура исследуемой системы
Здесь A(p) и B(p) – некоторые полиномы от p (свои для каждого варианта).
Характеристический полином системы равен знаменателю этой
передаточной функции
A ( p) =
(T1 p + 1)(T2 p + 1)(T3 p + 1). (2)
На основе приведенных выше критериев устойчивости можно заключить, что если коэффициенты T1, T2, T3 положительны, то система будет устойчивой. Анализ устойчивости замкнутой системы
можно осуществить с помощью критерия Найквиста.
1.2. Анализ управляемости и наблюдаемости
Понятия управляемости и наблюдаемости широко используются в современной теории автоматического управления. Линейная
система называется управляемой, если с помощью входного сигнала ее можно перевести из начала координат в любую точку пространства состояний. Система называется наблюдаемой, если по измерениям входного и выходного сигналов можно однозначно определить ее начальное состояние.
Анализ управляемости и наблюдаемости выполняется с помощью матриц управляемости и наблюдаемости или с помощью
51
грамианов управляемости и наблюдаемости. Те и другие строятся
по описанию в пространстве состояний
 = AX + bu, y = cX. X
(3)
Чтобы получить такое описание, воспользуемся структурной
схемой, приведенной на рис. 2. Выпишем операторные уравнения
для каждого блока схемы
=
x1
k1
k2
k3
=
u, x2
=
u, x3
(ax1 + x2 )
T1 p + 1
T2 p + 1
T3 p + 1
и преобразуем их к виду
T1 px1 = - x1 + k1u, T2 px2 = -x2 + k2u, T3 px3 = ak3 x1 + k3 x2 - x3 .
Выполняя обратное преобразование Лапласа, получаем систему
линейных дифференциальных уравнений первого порядка
T1x1 = - x1 + k1u,
T2 x2 = -x2 + k2u,
T3 x 3 = ak3 x1 + k3 x2 - x3 .
Кроме того, имеем алгебраическое уравнение для выходного сигнала y = x1–x3.
Переписывая эти уравнения в матричной форме (3), получаем
следующие выражения для матриц A, b, c:
0
0 
 a11
 b1 


 
A  0 a22
=
0=
b2  , c c1 0 c3  ,
 , b =
a31 a32 a33 
 0 
где ненулевые элементы матриц имеют вид
a11 =
-1 / T1, a22 =
- 1 / T2 , a31 =
ak3 / T3 , a32 =
k3 / T3 ,
a33 =
-1 / T3 , b1 =
k1 / T1, b2 =
k2 / T2 , c1 =
1, c3 =
-1.
Теперь можно перейти к анализу управляемости и наблюдаемости. Сформируем на основе матриц A, b, c две вспомогательные матрицы
R =  b, Ab, ..., An-1b  ,


52
 c 
 cA 
.
D=
 ... 
 n-1 
cA

Матрицы R и D называются соответственно матрицей управляемости и матрицей наблюдаемости системы. В пакете MATLAB их
можно построить с помощью команд ctrb и obsv.
Критерий управляемости. Для того чтобы система (3) была
управляемой, необходимо и достаточно, чтобы матрица управляемости имела полный ранг rankR = n.
Критерий наблюдаемости. Для того чтобы система (3) была наблюдаемой, необходимо и достаточно, чтобы матрица наблюдаемости имела полный ранг rankD = n.
В случае систем с одним входом и одним выходом матрицы R
и D квадратные, поэтому для проверки управляемости и наблюдаемости достаточно вычислить определители матриц R и D. Если они
не равны нулю, то матрицы имеют полный ранг.
Другой способ проверки управляемости и наблюдаемости опирается на вычисление грамианов управляемости и наблюдаемости. Так называются симметричные квадратные матрицы Wc и Wo,
определяемые равенствами
∞
=
Wc
T
eAt bbT eA t dt, Wo
∫=
0
∞
∫e
AT t T
c ceAt dt .
0
В пакете MATLAB их можно найти с помощью команд типа gram(sys, ‘c’), gram(sys, ‘o’). Необходимые и достаточные условия
управляемости и наблюдаемости имеют вид det Wc ≠ 0, det Wo ≠ 0.
1.3. Анализ минимальности моделей
Одной и той же передаточной функции Q(p) можно сопоставить целый класс эквивалентных реализаций в пространстве состояний, характеризуемых различными тройками матриц (A, b, c) разных, вообще
говоря, размеров. Реализация называется минимальной, если размер
ее матрицы A наименьший среди всех эквивалентных реализаций.
Поиск такой реализации имеет практический смысл, так как ее компьютерное моделирование требует меньших вычислительных затрат.
Для анализа минимальности реализации нужно проверить ее
управляемость и наблюдаемость.
Критерий минимальности. Для того чтобы реализация (3) была
минимальной, необходимо и достаточно, чтобы она была управляемой и наблюдаемой одновременно.
Таким образом, анализ минимальности конкретной реализации
сводится к проверке пары критериев rank R = n, rank D = n.
53
Если хотя бы один из рангов меньше n, то реализация неминимальна. Размерность эквивалентной минимальной реализации n0
определяется по формуле n0 = rank(RD).
Анализ минимальности с помощью грамианов управляемости и
наблюдаемости проводится аналогично. Критерием минимальности служит выполнение условия det(WcWo) ≠ 0, эквивалентного паре условий det Wc ≠ 0, det Wo ≠ 0.
Если исходное описание реализации оказалось неминимальным, то возникает задача перехода к минимальной реализации.
Чтобы решить ее, сначала перейдем от описания в пространстве состояний (3) к передаточной функции
Q( p) c=
( pE - A)-1 b
=
B( p)
.
A( p)
(4)
Далее нужно выделить общий множитель в числителе и знаменателе передаточной функции Q(p) и сократить на него. Эта процедура известна, как сокращение совпадающих нулей и полюсов
системы.
Отметим ряд соотношений между элементами матриц А, b, c
и коэффициентами передаточной функции, вытекающих из формулы (4). Знаменатель передаточной функции совпадает с характеристическим полиномом матрицы А
А(p) = det(pE – A) = pn + an–1pn–1 + ... + a1p + a0.
Его корни l1, ..., ln равны собственным числам матрицы А, а коэффициенты а0 и an–1 с точностью до знака равны ее определителю
и следу
a0 = det(–A ) = (–1)n l1 l2 ...ln,
–an–1 = trA = a11 + a22 + ... + ann = l1 + l2 + ... + ln.
(5)
Старший и младший коэффициенты числителя передаточной
функции
В(p) = bn–1pn–1 + ... + b1p + b0
удовлетворяют соотношениям
bn–1 = cb, b0/ a0 = –c A–1 b.
(6)
Величина k0 = b0/ a0 называется статическим коэффициентом
усиления. Она равна установившемуся значению переходной функции системы. В пакете MATLAB для вычисления статического коэффициента усиления имеется команда dcgain (от direct current
54
gain – коэффициент усиления по постоянному току). В ее основу положена формула k0 = Q(0).
Приведенные соотношения удобно использовать для контроля
вычислений.
Получение минимальной реализации в пакете MATLAB осуществляется командой minreal, для вычисления нулей и полюсов
можно использовать функции zero, pole, pzmap, zpk. Аргументом
во всех случаях служит исследуемая система sys, предварительно
сформированная командами ss или tf.
2. Задание по работе и содержание отчета
1. Нарисовать схему рис.2 для своего варианта и найти ее передаточную функцию по формуле (1).
2. Проверить устойчивость системы, используя алгебраический
и корневой критерии.
3. Получить для своего варианта схемы описание в пространстве
состояний вида (3) и выписать матрицы A, В, С.
4. Найти матрицы управляемости и наблюдаемости системы,
определить их ранги и сделать вывод об управляемости, наблюдаемости и минимальности системы.
5. Найти передаточную функцию по формуле (4), и сравнить ее
с полученной в п.1. Проверить выполнение соотношений (5), (6).
Тремя способами (по структурной схеме, передаточной функции и
описанию в пространстве состояний) найти статический коэффициент усиления системы.
6. Определить порядок минимальной реализации и найти ее передаточную функцию, выполнив сокращение нулей и полюсов.
Найти реакцию минимальной реализации на единичный скачок
и построить ее график.
7. Привести программы на языке MATLAB для выполнения пунктов 2–6 и краткое описание назначения и синтаксиса команд rank,
ss, tf, zpk, zero, pole, pzmap, ctrb, obsv, minreal, dcgain.
3. Порядок выполнения работы
1. В диалоговом режиме пакета MATLAB ввести матрицы A, b, c
и сформировать ss-описание системы sys = ss(A, b, c, 0). Используя
команды pole, eig, pzmap, найти полюсы системы и получить график их расположения на комплексной плоскости. Сделать вывод об
устойчивости. Построить диаграмму Найквиста системы и сделать
55
заключение об устойчивости системы, охваченной обратной связью
(рис. 1).
2. Найти матрицы управляемости и наблюдаемости, вычислить
их определители и ранги. Сделать вывод об управляемости, наблюдаемости и минимальности. Найти грамианы управляемости и наблюдаемости.
3. С помощью команд tf и zpk перейти к передаточной функции.
Двумя способами получить минимальную реализацию (сокращая
нули и полюса, и с помощью команды minreal). Сравнить переходные функции исходной и минимальной реализаций, а также их статический коэффициент усиления.
4. Набрать схемы исходной и минимальной реализаций в SIMULINK и сравнить их реакции на одинаковые входные сигналы.
4. Контрольные вопросы
1. Найти передаточные функции систем, заданных в пространстве состояний тройкой матриц:
 -1 0  1 
а) A = 
 ,   , c = 1 1 ,
 1 -2 0 
 -3 1 
1
б) A 
=
=
, b =

 , c 0 1 ,
 1 -2
1
 -4 1  1 
в) A = 
 ,   , c = 1 2 ,
 2 -3 0 
 -4 2  1
г) A = 
 ,   , c = 0 1 .
 1 -3 2
2. Найти ранги матриц управляемости и наблюдаемости для систем из п.1.
3. Нарисовать структурные схемы систем, матрицы которых
приведены ниже.
56
№
1
2
3
4
5
6
7
A
01
00
01
10
01
–1 0
01
11
10
02
10
01
11
11
bT
01
01
01
01
11
01
01
c
01
01
01
01
11
01
11
1
p+1
u
–
y
1
p+1
k
p+1
Рис. 3. Система управления третьего порядка
Определить, какие из них являются:
а) устойчивыми;
б) управляемыми;
в) наблюдаемыми;
г) минимальными.
4. Система управления задана структурной схемой, показанной
на рис. 3.
Требуется:
а) найти статический коэффициент усиления схемы;
б) найти передаточную функцию схемы и проанализировать ее
устойчивость;
в) найти описание схемы в пространстве состояний; построить
матрицы управляемости и наблюдаемости, сделать вывод о минимальности;
г) выяснить, при каких значениях коэффициента k схема будет
устойчивой, управляемой, наблюдаемой.
5. Варианты заданий
№ 1
2
3
4
5
6
7
8
9 10 11 12 13 14 15 16 17 18 19 20 21
k1 1
3
2
1
4
3
3
3
2
3
2
2
4
2
2
2
3
4
3
3
1
T1 1
4
2
3
1
3
3
2
2
4
1
1
2
3
5
2
2
4
1
1
4
k2 5
3
6
1
6
4
1
2
6
4
3
7
5
4
6
5
1
3
2
3
4
T2 1
4
2
3
1
3
3
2
2
4
2
2
2
4
5
2
2
4
1
1
4
k3 5
2
2
4
3
3
6
6
2
5
6
4
7
2
2
2
4
4
5
5
4
T3 4
2
1
2
4
5
3
3
4
3
5
5
2
2
5
2
4
2
2
4
3
a –2 –1 –1 –2 –1 –1 –1 –2 –1 –2 –2 –1 –2 –1 –1 –1 –2 –2 –2 –1 –1
57
Лабораторная работа №7
МОДЕЛИРОВАНИЕ РАЗНОСТНЫХ УРАВНЕНИЙ
Цель работы: ознакомиться с методикой моделирования разностных уравнений с помощью пакета MATLAB и SIMULINK.
1. Теоретические сведения [3]
Линейное однородное разностное уравнение второго порядка с постоянными коэффициентами имеет вид:
x(t + 2) + a1x(t + 1) + a0 x(t) =
f (t). (1)
Чтобы подчеркнуть дискретный характер изменения времени,
это уравнение часто записывают в форме
xk+2 + a1xk+1 + a0 xk =
fk .
Разностные уравнения третьего и более высоких порядков записываются аналогично.
Известны два основных метода решения линейных разностных
уравнений – с помощью характеристического полинома и с использованием z-преобразования, аналогичного преобразованию Лапласа.
1.1. Решение разностных уравнений
с помощью характеристического полинома
Будем искать решение однородного разностного уравнения в виде x(t) = zt, где z – некоторое число. Подставляя x(t) в разностное
уравнение (1) при f(t) = 0 и сокращая на zt, получаем характеристическое уравнение
z2 + a1z + a0 =
0.
Если его корни z1, z2 вещественные и различные, то общее решение имеет вид:
x=
(t) c1z1t + c2 z2t .
Если z1 = z2, то в решении появляется линейный множитель
x1 (=
t) (c1 + c2t)z1t .
58
В случае пары комплексно-сопряженных корней z1,2 = a ± ib решение может быть записано в вещественной форме
b
x(t) =rt (c1 sin jt + c2 cos jt), r= a2 + b2 , j= arctg ,
a
где r – модуль комплексного числа z1, а j – его аргумент.
Формулы для уравнений более высоких порядков выглядят так
же, просто увеличивается число слагаемых в решении.
Пример 1. Решим разностное уравнение
x(t + 2) – 5x(t + 1) + 6x(t) = 0.
Характеристическое уравнение имеет вид z2 –5z + 6 = 0.
Его корни вещественные и различные: z1 = 3, z2 = 2.
Общее решение: x(t) = c13t + c2 2t.
Пример 2. Решим разностное уравнение
x(t + 2) + 2x(t + 1) + 4x(t) = 0.
Характеристическое уравнение: z2 + 2z + 4 = 0. У него комплексные корни: z1,2 =-1 ± 3 i . Их положение на комплексной плоскости z = α + ib показано на рис. 1.
Модуль и аргумент корней можно найти непосредственно на
2p
рис. 1: r= 2, j=
.
3
2p
2p
Общее решение:
=
x(t) 2t (c1 cos t + c2 sin t).
3
3
Произвольные постоянные ci находят, задавая начальные условия. Пусть, например, в примере 2 заданы начальные условия
x(0) = 2; x(1) = –2. Записываем общее решение для t = 0 и t = 1:
2= 20 c1, -2= 21 (-c1 ⋅ 1 / 2 + c2 ⋅ 3 / 2).
2
z1
β
3
ρ
1
ϕ= 120°
–2
–1
ρ
z2
1
2 α
–1
– 3
–2
Рис. 1. Модуль и аргумент корней
59
Отсюда находим с1 = 2, с2 = 0. Следовательно, решение имеет вид
2p
x(t) = 2t +1 cos t.
3
Общее решение неоднородного разностного уравнения представляет собой сумму общего решения соответствующего однородного
уравнения и любого частного решения неоднородного уравнения.
Частное решение ищут в том же виде, что и правая часть, т. е. функция f(t) в уравнении (1):
– если f(t) – постоянная, то в виде константы;
– если f(t) – экспонента, то в виде экспоненты с тем же показателем;
– если f(t) = sin kt или cos kt, то в виде d1sinkt + d2coskt.
Коэффициенты d1 и d2 находят, подставляя частные решения
в разностное уравнение и приравнивая одноименные функции справа и слева.
Пример 3. Дано неоднородное разностное уравнение второго порядка
xk+2 + xk+1 + =
xk 6;
=
x0 1;
=
x1 1.
Находим корни характеристического полинома
z2 + z + 1= 0; z1,2=
-1 ± 3 i
p
; r= 1; j= .
2
3
Частное решение ищем в виде xчаст = d. Подставляя его в исходное уравнение, находим, что xчаст = 2.
Общее решение неоднородного уравнения получаем как сумму
частного решения и общего решения соответствующего однородного уравнения
p
p
xk =
2 + c1 sin k + c2 cos k.
3
3
Коэффициенты с1, с2 находим из уравнений
1=
2 + c1 3 / 2 + c2 / 2, откуда c1 =
- 3, c2 =
1.
1 = 2 + с2,
1.2. Решение разностных уравнений
с помощью z-преобразования
При описании дискретных систем и решении разностных уравнений широко применяется аппарат z-преобразования – это дискретный аналог преобразования Лапласа. Например, умножение
изображения F(p) на оператор Лапласа p соответствует дифферен60
цированию непрерывной функции f(t). Умножение изображения
F(z) на оператор z соответствует сдвигу функции f(t) (которая может
быть непрерывной, дискретной или решетчатой) на один такт.
Таким образом, если операторы р и 1/р – это операторы дифференцирования и интегрирования, то операторы z и z–1 – это операторы сдвига влево и вправо. С инженерной точки зрения оператор z–1
представляет собой элемент задержки.
Существуют также определенные параллели между изображениями функций F(p) и F(z). Например, изображению по Лапласу
F(p) = 1 соответствует дельта-функция f(t) = d(t), а изображению
F(z) = 1, соответствует единичный импульс. В том и другом случае
оригиналом является элементарное импульсное воздействие.
Краткая таблица z-преобразований других функций приведена
в Приложении.
Пусть дано разностное уравнение n-го порядка
y(t + n) + an–1y(t + n–1) + ... + a0y(t) = f(t)
(2)
с начальными условиями y(0) = y0; y(1) = y1; ...; y(n–1) = yn–1.
Алгоритм его решения с помощью z-преобразования содержит 4
шага.
Шаг 1. Применить z-преобразование к уравнению (2), заменяя
f(t) на F(z), y(t) на Y(z); y(t + 1) на z(Y(z)–y0) и т. д.
Шаг 2. Из полученного алгебраического уравнения выразить Y(z).
Шаг 3. Выполнить разложение Y(z) на простые дроби.
Шаг 4. Пользуясь таблицей, выполнить обратное z-преобразование.
Результатом будет искомое решение разностного уравнения.
Пример 4. Требуется решить разностное уравнение второго порядка
yn+2 – 5yn+1 + 6yn = un
с нулевыми начальными значениями y0, y1 и входным сигналом
un = 1.
Шаг 1. Применяем к нему z-преобразование
z2 Y (z) - 5zY (z) + 6 Y (z) =
U (z).
Шаг 2. Выражаем Y(z) и подставляем U (z) =
Y (z)
=
z
:
z -1
U (z)
1
z
=
⋅
.
2
2
z
z - 5z + 6 z - 5z + 6 - 1
61
Шаг 3. Представляем правую часть в виде суммы простых дробей с переменной z в числителе:
Y (z) =
z
0,5z
z
0,5z
.
=
+
(z - 1)(z - 2)(z - 3) z - 1 z - 2 z - 3
В отчете разложение делается методом неопределенных коэффициентов. В пакете MATLAB это можно сделать с помощью команды
residue.
Шаг 4. С помощью таблицы z-преобразований или команды
iztrans тулбокса SYMBOLIC пакета MATLAB находим оригиналы
каждого из слагаемых и складываем их:
yn = 0,5 – 2n + 0,5 · 3n.
1.3. Система линейных разностных уравнений
Матричная запись системы линейных однородных разностных
уравнений имеет вид
X(t + 1) = AX(t), X(0) = X0,
(3)
где X ∈ R n – вектор переменных, А – квадратная матрица постоянных коэффициентов, X0 – вектор начальных условий.
Решение этой системы может быть записано в компактной степенной форме
X(t) = At X0.
(4)
При работе в пакете MATLAB этот способ удобнее всего.
Пример 5. Найдем решение системы разностных уравнений второго порядка
x(t + 1) = 2x(t),
y(t + 1) = –2x(t) + 2y(t),
если=
x(0) 1=
; y(0) 20.
Матричное описание задачи имеет вид
 2 0
 x(0)   1 
X(t + 1) = AX(t), A = 
X(0) =
, =
  .
2
2


 y(0)  20 
Для получения решения воспользуемся степенной формулой.
Последовательно возводя матрицу А в степени 2, 3, ..., n (проделайте это!), находим
0
 1 0  1 
n  1
An = 2 
, X(t) = At X0 = 2t 

.

 -t 1  20 - t 
 -n 1 
62
Переходя к скалярной форме записи, получаем окончательный
ответ x = 2t, y = (20–t)2t.
Другой путь решения этой задачи связан с переходом от системы
двух разностных уравнений первого порядка к одному уравнению
второго порядка относительно y:
y(t + 2) – 4y(t + 1) + 4y(t) = 0.
Оно получается, если выразить из второго уравнения x(t), x(t + 1)
и подставить в первое.
1.4. Моделирование дискретных систем
в MATLAB и SIMULINK
В MATLAB при моделировании линейных систем с дискретным временем используются разностные уравнения, модели в пространстве состояний и дискретные передаточные функции. В ядре
MATLAB имеется команда filter, позволяющая рассчитать выходной сигнал фильтра по известному входному сигналу, заданному
массивом своих значений. Ее входными аргументами являются векторы а и b коэффициентов разностного уравнения, а также массив
значений сигнала и: y=filter(b,a,u). С помощью четвертого входного
аргумента можно задавать начальные условия фильтра.
Дискретная система в пространстве состояний, как и непрерывная система, задается четверкой матриц:
Xk+1 =
AXk + Buk , yk =
CXk + Duk , (5)
где Xk ∈ R n – вектор состояния; uk, yk – входной и выходной сигналы, A, B, C, D – постоянные матрицы.
Дискретной передаточной функцией линейной системы (5) называется отношение z-преобразования выходного сигнала к z-преобразованию входного сигнала при нулевых начальных условиях
Q(z) = Y(z)/U(z). Для систем с одним входом и одним выходом она
представляет собой отношение двух полиномов
Q(z) =
bn -1zn -1 + + b1z + b0
zn + an -1zn -1 + + a1z + a0
.
Дискретная передаточная функция системы, заданной описанием в пространстве состояний (5), может быть найдена по формуле
Q(z) = C(zE-A)–1B + D.
63
Дискретные модели в MATLAB задаются при помощи тех же
конструкторов ss, tf и zpk, что и непрерывные с той разницей,
что последним параметром задается частота дискретизации ts
(sampling time).
Пример 6. Рассмотрим дискретную систему, описываемую уравнениями
 -1 -1 / 2
1 
X
=
Xk +   uk ,=
yk 0 1 X,
k +1 

-1 
1 / 2
0 
при шаге дискретизации, равном 1.
Создадим соответствующую ss-модель с дискретным временем:
>> a=[-1 -1/2 ;1/2 -1]; b=[1;0]; c=[0 1]; d=0; ts=1.0; s1=ss(a,b,c,d, ts).
На дисплей будут выведены матрицы:
a=-1
0.5
-0.5
-1
b=1
0
c=0 1
d=0
и сообщение: Sampling time: 1.0 Discrete-time model.
Аналогично создаются tf и zpk-модели с дискретным временем.
Обращение к полям дискретной модели осуществляется так же, как
и у непрерывной, например для вывода матрицы c надо набрать код
s1.c. Время дискретизации хранится в поле ts.
Для моделирования дискретных систем используются те же команды (step, impulse, initial, lsim), что и в непрерывном случае. При
получении весовой функции дискретных систем вместо реакции на
дельта-функцию рассматривают реакцию на единичный импульс
вида
. Для этого используется та же команда impulse. Реакцию на единичный скачок (переходную функцию) получают при помощи команды step.
В тулбоксе SYMBOLIC имеются команды ztrans и iztrans для выполнения прямого и обратного z–преобразований.
Пример 7. Дискретная система задана уравнениями
yk+1 = xk , xk+1 = 5xk - 6yk + 7uk .
Найдем двумя способами (символьным и численным) передаточную функцию от входа u до выхода у.
В данном случае имеем следующие матрицы описания в пространстве состояний
y
 0 1
0 
=
X =
, A 
=
, B =


 , C [0 1].
x 
 -6 5
7 
64
1
(z-1)
z+0.5
z(z-0.5)
Discrete
Transfer Fcn
Discrete
Zero-Pole
1
y(n)=Cx(n)+Du(n)
x(n+1)=Ax(n)+Bu(n)
1+0.5z-1
Discrete Filter
Discrete State-Space
T
1
z-1
z
Discrete-Time Unit Delay
Integrator
Рис. 2. Дискретные блоки SIMULINK
Символьный способ. Дискретную передаточную функцию получаем по формуле Q = C(zE-A)–1B:
>>syms z;
>>А=[0 1;-6 5];В=[0; 7];С=[0 1];
>>Q=С*inv(z*eye(2)-А))*В Q=7*z/(z^2-5*z+6),
%ввод матриц
%передаточная функция
%результат
Численный способ. Используем команды тулбокса CONTROL:
>> A=[0 1;-6 5];B=[0; 7];C=[0 1];
>>sd=ss(A,B,C,0,1);
% Discrete-time model Sampling time: 1,.
>> s=tf(sd)
Transfer function:
7z
---------------z^2 - 5 +6
В обоих случаях получаем одинаковый результат:
Q(z) =
7z
2
z - 5z + 6
.
Моделирование дискретных систем в SIMULINK производится
так же, как и непрерывных. На рис. 2 приведены основные блоки из
библиотеки Discrete.
Блок задержки Unit delay позволяет строить схемы моделирования методом Кельвина. Блоки Discrete Transfer Fcn (дискретная передаточная функция) и Discrete Filter (дискретный фильтр) определяют
рациональную функцию отношения полиномов от оператора z и обратного оператора 1/ z соответственно.
2. Задание по работе и содержание отчета
Работа состоит из двух частей.
1. В первой части исследуется дискретная система, заданная разностным уравнением yn + 2 + byn + 1 + ayn = un с нулевыми начальными условиями y0 = 0; y1 = 0.
65
Требуется для своего варианта найти тремя способами ее реакцию на входной сигнал un = 1:
а) последовательно рассчитав точки y2, …, y5;
б) решив разностное уравнение;
в) используя z-преобразование.
Привести все числовые выкладки по решению разностного уравнения, дискретную передаточную функцию, ее разложение на
простые дроби, график yn, схему моделирования разностного уравнения в SIMULINK, построенную на элементах задержки 1/z и программы для вычисления в MATLAB.
2. Во второй части работы рассматривается популяция (рыб, животных), разбитых на три возрастные группы – младшего, среднего
и старшего возрастов. Заданы величины р1, р2 – вероятности дожития особями каждой возрастной группы до следующего возраста,
и числа а1, а2, а3, характеризующие среднюю плодовитость каждой
возрастной группы.
Нужно выяснить, как изменяется со временем численность возрастных групп, и каково будет их процентное соотношение через достаточно большое время.
Возрастная структура популяции описывается матричным уравнением
x 
 a1 a2 a3 
 


X(=
k + 1) AX(k=
), X  y =
, A  p1 0 0  , (6)
 z 
 0 p2 0 
где x, y, z – численности трех возрастных групп.
Требуется рассчитать численности возрастных групп для n = 1, 2, 3
при начальных условиях x0 = y0 = z0 = 10. Теоретически найти процентное соотношение численностей x, y, z для больших n – оно определяется
соотношением компонент главного собственного вектора матрицы А.
Отчет должен содержать все необходимые расчеты, программы
для вычисления в MATLAB и схемы моделирования в SIMULINK.
Для определения собственных векторов матрицы А можно использовать символьный вариант команды eig.
3. Порядок выполнения
1. Получить решение разностного уравнения и построить графики:
а) с помощью команд ztrans, iztrans, ezplot тулбокса SYMBOLIC;
б) с помощью команд tf и step тулбокса CONTROL;
в) путем моделирования в SIMULINK;
66
2. Рассчитать абсолютные и относительные численности возрастных групп популяции:
а) с помощью степенной формулы (4);
б) с помощью команды initial.
Построить графики в обычном и логарифмическом масштабах.
4. Контрольные вопросы
Даны разностные уравнения:
1. Даны разностные уравнения:
1. x(t + 2) –5x(t + 1) + 6x(t) = 0,
2. x(t + 2) –3x(t + 1) + x(t) = 0,
3. x(t + 2) – 3x (t + 1) + 2x(t) = 0,
4. yk+1 + yk + 5yk–1 + 3yk–2 = 1,
5. yn+2 = yn+1 – yn, y1 = 1, y2 = 0,
6. xk+3 + 7xk+2 + 15xk+1 + 9xk = 0,
7. xk+2–3xk+1 + 2xk = 0,
8. x(t + 3) + 7x(t + 2) + 15x(t + 1) + 9x(t) = 0,
9. x(t + 4) + 2x(t + 2) + x(t) = 0,
10. yn+3 – 5yn+2 + 8yn+1 – 4yn = 0; y0 = 0, y1 = 2, y2 = 10,
11. yk+1 - yk + 2yk-1 =
0,
12. 4yk+1 + 4yk + yk-1 =
0,
13. yk+2 + y=
0
,
y
=
2
, y=
k
0
1 1.
14. yk+1 - 4yk + yk-1 + 6yk-2= 0, y0= 6, y1= 12, y2= 276.
Требуется для уравнения, указанного преподавателем:
а) Решить его с помощью характеристического уравнения;
б) Решить его при помощи z-преобразования;
в) Нарисовать схему моделирования и от нее перейти к матричной форме записи (3).
2. Решить систему разностных уравнений двумя способами: по
формуле (4) и путем перехода к одному уравнению
-3xê + y ê , x0 =
1,
xêí =
а) 
xê - 3y ê , y0 =
1.
yêí =
x(t) + 3y(t)
x(t + 1) =
б) 
y(t + 1)= x(t) + y(t).
=
yn +1 - x=
n 0, y
0 2,
в) 
=
=
n 0, x
0 0.
xn +1 - y
67
3. Решить систему Xn +1 = AXn , X0 для трех вариантов матрицы А:
 -1 2 -4 


а) A0 =
 -8 -3 2  ,
 -2 -4 6 
=
в) A2
T
A=
1 ; X0
0 2 1 


б) A1 = 1 0 2  ,
1 2 0 
20 
 
 0 .
 0 
5. Варианты заданий
Варианты коэффициентов a и b разностного уравнения yn+2 +
+ byn+1 + a yn = 1
№
1
2
3
4
5
6
7
8
9
10
b
–1,3
–0,7
0,2
–2
–1,6
1,6
–1,6
1,9
–2,8
2,5
a
0,3
–0,6
–0,48
0,96
0,6
0,55
0,48
0,6
1,8
1
№
11
12
13
17
18
19
20
b
0,5
2
–2,88
a
0,27
1,1
1,5
–0,72 1,28
0,4
0,7
14
15
16
–2
2,88
–0,5
0,72 –1,28
1,1
1,6
0,25
0,37
0,66
Варианты элементов матрицы А системы разностных уравнений (6) (а1 = 0)
№
1
2
3
4
5
6
7
8
а2
9
9
2
9 18 6
3
6 12
9 10 11 12 13 14 15 16 17 18 19 20
27 27 27
27 12 6
2 4 2
2 24
27 27 27 27
4 8 4 2
а3 24 12 6 12 16 8
4
8 36 18 9 18 32 16 8 16 18 9 4,5 9
р1 1
3
2
3
1
3
1
6
2
9
4
9
2
9
1
9
1
2
1
1
2
1
4
4
9
8
9
4
9
2
9
1
4
1
8
р2 1
4
1
4
1
1
1
6
1
6
2
3
2
3
3
8
3
8
3
2
3
2
1
3
1
3
4
3
4 3 3 3
3 16 16 4
3
4
68
1
4
1
2
Библиографический список
1. Мироновский Л.А. Моделирование линейных систем: учебное
пособие. СПб.: ГУАП, 2009.
2. Мироновский Л. А., Петрова К. Ю. Введение в MATLAB: Учебное пособие. СПб.: ГУАП, 2006.
3. Мироновский Л. А. Моделирование разностных уравнений:
учебное пособие. СПб.: ГУАП, 2004.
69
ПРИЛОЖЕНИЕ
Некоторые математические формулы
2
1. Корни квадратного уравнения x2 + 2bx + c =
0, x1,2 =-b ± b - c.
2
2. a sin t + b cos=
t c sin ( t + j ), c=
a2 + b2 , tgj
= b / a.
c
b
ϕ
a
Производные от элементарных функций
1. c′ = 0, ( )′ = aeat , 4. eat
(
)
(
)
( )′ = 2t, ( )′ = ktk-1,
2. t2
3. tk
5. ( sin kt )′ = k cos kt, 6. ( cos kt )′ = -k sin kt,
′
7. eat sin kt =eat ( a sin kt + k cos kt ), ′
8. eat cos kt =eat ( a cos kt - k sin kt ).
Интегралы простейших функций
1. ∫ ad=
t at + c, k
3. ∫ t=
dt
t
tk+1
+ c, k +1
t
5. ∫ et d=
t et = et - 1, 0
0
2. ∫ td=
t
t2
+ c, 2
4. ∫ et d=
t et + c, dt
6. ∫ eat
=
1 at
e + c,
a
1
1
ktdt
sin kt + c,
7. ∫ sin ktdt =
=
- cos kt + c, 8. ∫ cos
k
k
9.
70
dt
∫=
t
ln t + c.
Таблица преобразований Лапласа
Оригинал f(t)
Дельтафункция
Изображение
F(р)
1
1(t)
t
t2
e–at
1
1
2
1
p
2
3
p
e–at sinωt
Оригинал f(t)
p
sinωt
ω
2
p +ω
te–at
ω
p+a
1
Изображение
2
2
2
2
F(р)
( p + a) + ω ( p + a) + ω ( p + a)2
p
2
p+a
e–at cosωt
cosωt
2
p + ω2
tsinωt
tcosωt
2 pω
p2 - ω2
( p2 + ω2 )2
( p2 + ω2 )2
В тулбоксе SYMBOLIC пакета MATLAB имеются команды
laplace и ilaplace для выполнения прямого и обратного преобразований Лапласа.
Например, набрав код >> syms t; laplace(t*exp(-t),’p’), получим
ans=1/(p+1)^2.
Таблица z-преобразований
Единичный
импульс
Изображение
F(z)
1
Изображение F(z)
Оригинал f(t)
Оригинал f(t)
sinωt
t2
t
1
z
z
z(z + 1)
z -1
(z - 1)2
(z - 1)
z sin ω
e–at
z
z
z-a
z(z - cos ω)
2
z - 2z cos ω + 1 z - 2z cos ω + 1 z - 2az cos ω + a
Оригинал f(t)
shwt
chwt
Изображение
F(z)
zshω
z(z - chω)
2
2
z - 2zchω + 1 z - 2zchω + 1
a
(z - a)2
atcoswt
az sin ω
2
t at
az
z-e
atsinwt
cosωt
2
3
at
z(z - a cos ω)
2
2
z - 2az cos ω + a2
f(t + 1)
f(t + 2)
z (F(z) – f0)
z2(F(z) –
– f0 – z–1f1)
71
В тулбоксе SYMBOLIC пакета MATLAB имеются команды ztrans
и iztrans для выполнения прямого и обратного z-преобразований.
Например, набрав код >>syms t; ztrans(t^2), получим
ans = z*(z + 1)/(z–1)^3.
Таблица функций e–t, sin t, cos t
t
0,0
0,1
0,2
0,3
0,4
0,5
0,6
0,7
0,8
0,9
e–t
1,000 0,905 0,819 0,741 0,670 0,606 0,549 0,497 0,449 0,407
sint
0,000 0,100 0,199 0,296 0,389 0,479 0,565 0,644 0,717 0,783
cost
1,000 0,995 0,980 0,955 0,921 0,878 0,825 0,765 0,697 0,622
t
1,0
1,6
1,7
1,8
0,202
0,183
0,165
sint 0,842 0,891 0,932 0,964 0,985 0,997 1,000 0,9996
0,9917
0,9738
e–t
1,1
1,2
1,3
1,4
1,5
1,57
0,368 0,333 0,301 0,272 0,247 0,223 0,208
cost 0,540 0,454 0,362 0,268 0,170 0,071 0,001 –0,0292 –0,1288 –0,2272
t
2,0
2,1
2,2
2,3
2,4
2,5
0,150
0,135
0,122
0,111
0,100
0,091
0,082
sint
0,9463
0,9093
0,8632
0,8085
0,7457
0,6755
0,5985
cost
–0,3233
–0,4161
–0,5048
–0,5885
t
2,6
2,7
2,8
2,9
3,0
3,1
3,2
e–t
0,051
0,050
0,045
0,041
0,050
0,045
0,041
sint
0,5155
0,4274
0,3350
0,2392
0,1411
0,0416
0,0584
cost
72
1,9
e–t
–0,8569 –0,9041
–0,9422
–0,6663 –0,7374
–0,8011
–0,9710 –0,9900 –0,9991 –0,9983
73
Анализ линейных систем
pole, zero, pzmap, bode, nyquist, ltiview,
ctrb, obsv, minreal
Моделирование
step, impulse, initial, lsim,
filter,
ode23, ode45
Линейные системы
ss, tf, zpk, ss2tf, ss2zp
Графика
Генераторы и осцилографы
const, step, sine wave,
scope, XY Graph
Линейные блоки
sum, gain, abs, integrator, transfer fcn,
state-space
Блоки SIMULINK
mux, in, out, to Workspace
Вспомогательные блоки
solve, dsolve, syms,
laplace, ilaplace,
ztrans, iztrans
Символьные операции
plot,, subplot, ezplot,
zeros, ones, eye, diag det, rank,
mesh, surf, ezsurf, grid,
inv, eig, svd
hold, figure
roots, poly, polyval, conv,
residue
abs, angle, sqrt, exp, sin,
cos, tan, sinh, cosh, stepfun
Матрицы
Полиномы
Функции
Команды MATLAB
Некоторые команды MATLAB и блоки SIMULINK
СОДЕРЖАНИЕ
Предисловие............................................................... Лабораторная работа № 1.
Модели линейных блоков.............................................. Лабораторная работа № 2.
Моделирование дифференциальных уравнений................ Лабораторная работа № 3.
Моделирование следящей системы................................. Лабораторная работа № 4.
Моделирование двойного маятника................................ Лабораторная работа № 5.
Моделирование колебаний струны................................. Лабораторная работа № 6.
Анализ системы управления.......................................... Лабораторная работа №7.
Моделирование разностных уравнений........................... Приложение................................................................ 74
3
5
15
22
29
38
49
58
70
Для заметок
75
Учебное издание
Мироновский Леонид Алексеевич
МОДЕЛИРОВАНИЕ ЛИНЕЙНЫХ СИСТЕМ
в MATLAB и SIMULINK
Лабораторный практикум
Публикуется в авторской редакции
Компьютерная верстка Н. Н. Караваевой
Сдано в набор 16.11.16. Подписано к печати 09.01.17.
Формат 60×84 1/16. Бумага офсетная. Усл. печ. л. 4,42.
Уч.-изд. л. 4,63. Тираж 50 экз. Заказ № 8.
Редакционно-издательский центр ГУАП
190000, Санкт-Петербург, Б. Морская ул., 67
76
Документ
Категория
Без категории
Просмотров
0
Размер файла
2 111 Кб
Теги
0bbf6341ac, mironovskiy
1/--страниц
Пожаловаться на содержимое документа