close

Вход

Забыли?

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

?

Mironovskii

код для вставкиСкачать
Министерство образования и науки российской федерации
Федеральное государственное автономное образовательное
учреждение высшего профессионального образования
САНКТ-ПЕТЕРБУРГСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ
АЭРОКОСМИЧЕСКОГО ПРИБОРОСТРОЕНИЯ
ТЕОРИЯ ОПТИМАЛЬНОГО
УПРАВЛЕНИЯ
Методические указания
к выполнению лабораторных работ
А1
u
v
Санкт-Петербург
2012
А2
Составитель Л. А. Мироновский
Рецензенты: кафедра информационных систем СПбГУАП;
кандидат технических наук Г. С. Бритов
Содержатся указания к выполнению лабораторных работ по методам оптимизации на персональных ЭВМ с использованием программного пакета MATLAB.
Методические указания предназначены для проведения лабораторных работ по курсу «Теория оптимального управления», со студентами дневного обучения по направлению 230100.68 «Информатика
и вычислительная техника».
Подготовлены к изданию кафедрой вычислительных систем и сетей и рекомендованы редакционно-издательским советом СанктПетербургского государственного университета аэрокосмического
приборостроения.
Редактор А. В. Подчепаева
Компьютерная верстка Н. Н. Караваевой
Сдано в набор 13.09.12. Подписано к печати 12.10.12. Формат 60×84 1/16.
Уч.-изд. л. 2, 56. Усл. печ. л. 2,38. Бумага офсетная. Тираж 100 экз. Заказ № 532.
Редакционно-издательский центр ГУАП
190000, Санкт-Петербург, Б. Морская ул., 67
© Санкт-Петербургский государственный
университет аэрокосмического
приборостроения (ГУАП), 2012
Лабораторная работа № 1
СИНТЕЗ ПИД-РЕГУЛЯТОРА
Цель работы: освоение методики настройки ПИД-регулятора
и исследование возможностей его применения для оптимизации переходной характеристики системы в пакете MATLAB.
1. ТЕОРЕТИЧЕСКИЕ СВЕДЕНИЯ
1.1. Структура системы управления
Важной задачей теории управления является синтез регуляторов, обеспечивающих заданные динамические характеристики
конструируемой системы.
В работе исследуется так называемый ПИД-регулятор, в состав
которого входят пропорциональное (П), интегрирующее (И) и дифференцирующее (Д) звенья для управления объектом. Рассмотрим
систему с единичной обратной связью, изображенную на рис. 1.
ПИД-регулятор представляет собой параллельное соединение
указанных звеньев (рис. 2), поэтому его передаточная функция выглядит следующим образом:
KP +
KI
K s2 + KP s + KI
+ KD s = D
,
s
s
(1)
где KP – пропорциональный коэффициент усиления; KI – интегральный коэффициент усиления; KD – дифференциальный коэффициент усиления.
хx
e
+
Регулятор
u
–
Объект
управления
·
y
Рис. 1
3
KI/s
e
Kp
u
K0/s
Рис. 2
Работа замкнутой системы, показанной на рис. 1, происходит
следующим образом. Переменная e представляет ошибку слежения как разницу между задаваемым входным значением х и текущим выходом y. Этот сигнал ошибки поступает в ПИД-регулятор,
который вычисляет производную
и интеграл ошибки. Сигнал u после регулятора имеет вид
u = KP e + KI ò edt + KD
de
.
dt
(2)
Он поступает на вход объекта и приводит к изменению его выхода y. Это значение выхода вновь поступает на вход для нахождения
нового сигнала ошибки. Регулятор получает этот сигнал и снова интегрирует и дифференцирует его. Этот процесс непрерывно повторяется.
1.2. Влияние коэффициентов ПИД-регуляторов
на характеристики переходного процесса
Увеличение коэффициента усиления KP пропорционального регулятора сокращает время нарастания выходного сигнала и уменьшает, но не сводит к нулю, установившуюся ошибку. Интегрирующий регулятор полностью устраняет установившуюся ошибку, но
сильно ухудшает переходную характеристику. Дифференциальный
регулятор увеличивает устойчивость системы, уменьшает перерегулирование и улучшает переходную характеристику.
Влияние каждого регулятора в замкнутой системе показано
в табл.1.
Таблица 1
Тип
звена
4
Время
Нарастания
Перерегулирование
Время переходного
процесса
Статическая
ошибка
KP
Уменьшает
Увеличивает
Слабо влияет
Уменьшает
KI
Уменьшает
Увеличивает
Увеличивает
Исключает
KD
Слабо влияет
Уменьшает
Уменьшает
Слабо влияет
Отметим, что эти зависимости
x
могут быть не очень точными, по.
тому что степень влияния каждоbx
го из коэффициентов KP, KI, и KD
F
зависит от значений других коM
k
эффициентов. Изменение одного
из них может изменить эффект
остальных двух. Поэтому табл. 1
может быть использована тольРис. 3
ко как рекомендация при выборе
величин KI, KP и KD.
Пример. Предположим, что у нас есть подвижная масса М, пружина жесткости k и демпфер с коэффициентом демпфирования b
(рис. 3).
Обозначим смещение массы от положения равновесия через х,
тогда величины х̇ и х̇˙ будут характеризовать ее скорость и ускорение. Уравнение моделирования этой системы при наличии управляющей силы F записывается на основе второго закона Ньютона:
 + bx + kx = F.
Mx
Применим к обеим частям этого уравнения преобразование Лапласа, полагая начальные условия нулевыми:
Ms2 X(s) + bsX(s) + kX(s) = F (s).
Передаточная функция от входа F(s) до смещения X(s) имеет вид
X(s)
1
=
.
2
F (s) Ms + bs + k
(3)
Положим M = 1 кг, b = 10 н с/м, k = 20 н/м, F(s) = 1 н и подставим
эти величины в передаточную функцию:
X(s)
1
=
.
2
F (s) s + 10s + 20
(4)
Попытаемся с помощью компьютерного моделирования выяснить, как нужно изменять коэффициенты KP, KI и KD, чтобы обеспечить:
– увеличение скорости переходного процесса;
– минимизацию перерегулирования;
– уменьшение установившейся ошибки.
5
1.3. Переходная характеристика разомкнутой системы
Получим реакцию разомкнутой системы на ступенчатое воздействие в пакете MATLAB. Создадим m-файл, в который запишем
следующие команды:
num=1; den=[1 10 20]; sys=tf(num, den); step(sys)
Запустив его в командном окне MATLAB, получим график, показанный на рис. 4.
Статический коэффициент усиления передаточной функции равен 1/20, поэтому установившееся значение выходной величины
при единичном ступенчатом воздействии равно 0,05. Это соответствует очень большой установившейся ошибке – 0,95. Кроме того,
постоянная времени разомкнутой системы около одной секунды,
а время переходного процесса около 1,5 секунд. Требуется разработать регулятор, который уменьшит постоянную времени, сократит
время переходного процесса и устранит статическую ошибку.
1.4. Пропорциональный регулятор
Из табл. 1 видно, что пропорциональный регулятор уменьшает
постоянную времени, но увеличивает перерегулирование. Передаточная функция замкнутой системы с пропорциональным регулятором имеет вид:
KP
X(s)
= 2
.
F (s) s + 10s + (20 + KP )
(5)
Open-Loop Step
0.05
Displacement (m)
0.04
0.03
0.02
0.01
0
0
0.5
1
Time (sec)
Рис. 4
6
1.5
2
Зададим пропорциональный коэффициент KР = 300 и изменим
m-файл следующим образом:
KР=300; num=[KР]; den=[1 10 20+KР]; t=0:0.01:2;
sys=tf(num,den); step(sys, t)
Запустив его в командном окне MATLAB, увидим следующий
график (рис. 5).
Полученный график показывает, что пропорциональное звено
уменьшает время переходного процесса и статическую ошибку, но
приводит к появлению перерегулирования.
Заметим, что для нахождения передаточной функции замкнутой системы по известной передаточной функции разомкнутой системы можно использовать команду MATLAB feedback. Построение того же графика с ее помощью выглядит следующим образом:
num=1; den=[1 10 20]; Kp=300;
sys=tf(num, den), sysF=feedback(Kp*sys, 1),
step(sysF)
1.5. Пропорционально-дифференциальный регулятор
Перейдем к рассмотрению ПД-управления. Из табл. 1 видно, что
коэффициент KD уменьшает величину и время перерегулирования.
Передаточная функция замкнутой системы с ПД-регулятором имеет вид
KD s + KP
X(s)
. = 2
(6)
F (s) s + (10 + KD )s + (20 + KP )
Closed-Loop Step: Kp=300
1.4
Displacement (m)
1.2
1
0.8
0.6
0.4
0.2
0
0
0.5
1
Time(sec)
1.5
2
Рис. 5
7
Зададим пропорциональный коэффициент KР = 300 и дифференциальный коэффициент KD = 10. Сформируем и выполним следующий m-файл:
KР=300; KD=10; num = [KD KP]; den = [1 10+KD 20+KP];
t = 0 :0.01 :2 ; step(num,den,t).
Результатом будет график переходной характеристики, показанный на рис. 6.
Из него видно, что дифференциальный регулятор уменьшает
ошибку отклонения и время переходного процесса и слабо влияет
на постоянную времени и статическую ошибку.
1.6. Пропорционально-интегральный регулятор
Прежде, чем перейти к ПИД-регулятору, рассмотрим ПИ-регулятор. Из табл. 1 следует, что интегрирующий регулятор устраняет статическую ошибку, но при этом увеличивает перерегулирование и время переходного процесса. Передаточная функция системы
с ПИ-регулятором имеет вид
KP s + KI
X(s)
.
=
F (s) s3 + 10s2 + (20 + KP )s + KI
Уменьшим коэффициент KP до 30, а KI установим равным 70.
Создадим новый m-файл, содержащий следующие команды:
KP=30; KI=70; num=[KP KI]; den=[1 10 20+KP KI];
t=0:0.01:2; step(num,den,t)
Closed-Loop Step: Kp = 300, Kd = 10
1.4
Displacement (m)
1.2
1
0.8
0.6
0.4
0.2
0
0
0.5
1
Time (sec)
Рис. 6
8
1.5
2
Запустив m-файл в командном окне MATLAB, получим график,
показанный на рис. 7.
Мы уменьшили пропорциональный коэффициент KР, поскольку интегрирующий регулятор так же уменьшает постоянную время
и увеличивает перерегулирование, как и пропорциональный регулятор (двойной эффект). Полученный график показывает, что интегральный регулятор полностью устранил статическую ошибку.
1.7. Пропорционально-интегральнодифференциальный регулятор
Теперь рассмотрим ПИД-регулятор. Передаточная функция замкнутой системы с ПИД-регулятором следующая:
KD s2 + KP s + KI
X(s)
.
= 3
F (s) s + (10 + KD )s2 + (20 + KP )s + KI
(7)
После нескольких шагов компьютерного моделирования с помощью метода проб и ошибок приходим к значениям KР = 350,
KI = 300 и KD = 50, обеспечивающим желаемую реакцию (рис. 8).
Для построения графика использовался следующий m-файл:
Kp=350; Ki=300; Kd=50;
num=[Kd Kp Ki]; den=[1 10+Kd 20+Kp Ki];
t=0:0.01:2; step(num,den,t)
Из графика видно, что мы получили систему без перерегулирования, с быстрым временем нарастания и без статической ошибки.
Closed-Loop Step: Kp = 30, Ki = 70
1.4
Displacement (m)
1.2
1
0.8
0.6
0.4
0.2
0
0
0.5
1
Time (sec)
1.5
2
Рис. 7
9
Closed-Loop Step: Kp = 350, Ki = 300, Kd = 5500
1
Displacement (m)
0.8
0.6
0.4
0.2
0
0
0.5
1
Time (sec)
1.5
2
Рис. 8
1.8. Основные приемы синтеза ПИД-регулятора
При разработке ПИД-регулятора для заданной системы можно
рекомендовать следующую процедуру его синтеза, обеспечивающую получение желаемого результата.
Шаг 1. Получите реакцию разомкнутой системы и определите,
что должно быть улучшено.
Шаг 2. Добавьте пропорциональное звено для уменьшения постоянной времени.
Шаг 3. Добавьте дифференцирующее звено для уменьшения перерегулирования.
Шаг 4. Добавьте интегрирующее звено для устранения статической ошибки.
Шаг 5. Варьируйте каждый из коэффициентов KР, KI и KD до
тех пор, пока не получите желаемый результат. При этом можно использовать приведенную ранее табл. 1, которая показывает, на что
влияют отдельные элементы регулятора.
Следует иметь в виду, что совсем не обязательно применять все
три части регулятора (пропорциональную, дифференциальную и
интегральную). Например, если ПИ-регулятор дает достаточно хорошую реакцию (как в рассмотренном примере), то не нужно добавлять дифференцирующее звено. Надо стараться создавать по возможности более простой регулятор.
10
2. ЗАДАНИЕ ПО РАБОТЕ И СОДЕРЖАНИЕ ОТЧЕТА
Требуется найти коэффициенты ПИД-регулятора для заданного
варианта системы, приведенной на рис. 3.
Для исследования это системы необходимо:
1) теоретически найти реакцию разомкнутой системы на ступенчатое входное воздействие, а также реакцию системы с единичной
обратной связью. Построить графики y(t);
2) получить аналогичные графики для случаев П-регулятора,
И-регулятора, Д-регулятора и их сочетаний;
3) привести MATLAB-программы и SIMULINK-схемы для моделирования разомкнутой и замкнутой систем.
Отчет должен содержать вывод формул (3)–(7), аналитические
выражения для переходной функции во всех случаях и построенные по ним графики (без использования компьютера!).
3. ПОРЯДОК ВЫПОЛНЕНИЯ
Работа выполняется сначала в диалоговом режиме в командном
окне MATLAB, а потом путем моделирования в SIMULINK. Результаты сравниваются.
В обоих случаях сначала следует получить переходную характеристику разомкнутой системы, затем наблюдать ее изменение при
введении пропорциональной, интегрирующей и дифференцирующей частей ПИД-регулятора. Получаемые результаты и графики
следует сравнивать с теоретическими, приведенными в отчете.
При настройке ПИД-регулятора в SIMULINK наблюдать реакцию системы на синусоидальный входной сигнал.
Примечание. В последних версиях SIMULINK ПИД-регулятор
введен как стандартный блок. Попробуйте использовать его при выполнении работы.
4. КОНТРОЛЬНЫЕ ВОПРОСЫ
1. Опишите процедуру настройки ПИД-регулятора.
2. Как меняются нули и полюсы передаточной функции при введении П-регулятора? И-регулятора? Д-регулятора?
3. Нарисуйте на комплексной плоскости траектории корней при
изменении коэффициента КР от 0 до 300.
4. Какие параметры или сочетания параметров не меняются при
введении П-, И-, Д-регуляторов?
5. В чем отличие ПИД-регулятора от модального регулятора? Что
у них общего?
11
5. ВАРИАНТЫ ЗАДАНИЙ
12
№
1
2
3
4
5
6
7
8
9
10
М
1
1
1
1
1
1
1
1
2
2
b
0,1
0,2
0,3
0,4
0,1
0,2
0,3
0,4
0,1
0,2
k
1
1
1
1
4
4
4
4
4
4
№
11
12
13
14
15
16
17
18
19
20
M
2
2
1
1
1
1
2
2
2
2
B
0,3
0,4
0,1
0,2
0,3
0,4
0,1
0,2
0,3
0,4
k
4
4
9
9
9
9
9
9
9
9
Лабораторная работа № 2
МОДАЛЬНОЕ УПРАВЛЕНИЕ
Цель работы: освоить методику расчета модального управления
с обратной связью в линейных динамических системах и исследовать возможности модального управления для изменения динамических характеристик системы на основе использования пакета
MATLAB.
1. ТЕОРЕТИЧЕСКИЕ СВЕДЕНИЯ
1.1. Принцип модального управления
Основной задачей теории управления является отыскание законов управления, обеспечивающих заданные динамические характеристики конструируемой системы. Чаще всего для решения этой
задачи используется принцип обратной связи, когда управляющий
сигнал формируется в зависимости от переменных состояния xi(t)
или от текущих значений выходного сигнала системы y(t). Для линейных динамических систем хорошие результаты получаются,
когда закон управления является линейной функцией состояния
или выходного сигнала.
Рассмотрим задачу синтеза линейного закона управления для
полностью управляемой системы, заданной описанием в пространстве состояний
 (t) = AX(t) + bu(t),
X
(1)
y
(
t
)
=
c
X
(
t
).
(2)
где X(t) Î R n – вектор состояния; u(t), y(t) Î R – входной и выходной
сигналы; A – (n × n)-матрица коэффициентов; b – (n × 1)-матрица
входов; с – (1 × n )-матрица выходов.
Предположим, что для использования в управлении доступны
все переменные состояния x1, ..., xn системы.
Обозначим через λ1,, λn собственные числа матрицы А (полюсы системы). Именно от них зависят такие динамические свойства системы, как устойчивость, быстродействие, колебательность
и другие. Цель модального управления состоит в том, чтобы изменить значения λ i , сделав их равными некоторым желаемым (назначаемым) числам µ1,, µn . Чтобы добиться этого, введем для
13
системы (1) линейное управление u = -k1x1 -  - kn xn + k0v или
в векторной форме записи
u(t) = -kX(t) + k0v(t), (3)
где k = [k1, ..., kn] – вектор-строка параметров обратной связи, v – новая входная переменная.
Подставляя это управление в уравнение (1), получаем
 (t) = AX(t) - bkX(t) + bk v(t) = (A - bk)X(t) + bk v(t).
X
0
0
(4)
Известно, что если система (1) полностью управляема, то выбором коэффициентов ki можно добиться любого желаемого расположения собственных чисел µi матрицы (A–bk) замкнутой системы,
независимо от расположения собственных чисел λi матрицы А разомкнутой системы. Тем самым могут быть решены задачи:
а) обеспечения устойчивости управляемой замкнутой системы,
если разомкнутая система была неустойчивой;
б) улучшения динамических характеристик замкнутой системы, если эти характеристики (например, величина перерегулирования, коэффициент затухания и др.) для разомкнутой системы не
удовлетворяли конструктора.
Поскольку собственные числа λi матрицы А непосредственно определяют собственные движения («моды») системы eλ it , линейное
управление, обеспечивающее заданный набор собственных чисел
замкнутой системы, получило название «модальное управление».
1.2. Методика расчета параметров модального управления
Итак, необходимо для линейной управляемой системы (1), матрица А которой имеет набор собственных чисел {λ i }, i = 1, n, , построить управление вида (3) такое, чтобы матрица (A – bk) замкнутой системы имела бы желаемый набор собственных чисел {µ i }, i = 1, n. .
Напомним, что собственные числа матриц A и A–bk представляют собой корни характеристических полиномов этих матриц, которые вычисляются по формулам
F = |pE–A| и Fзам = |pE–A+bk|.
Поскольку собственные числа матрицы однозначно определяют коэффициенты ее характеристического полинома, задача может быть сформулирована следующим образом: для
управляемой системы (1) с характеристическим полиномом
F = pn + αn-1 pn-1 + ... + α 0 найти вектор k коэффициентов обратной связи такой, чтобы замкнутая система (4) имела бы характери14
n
n-1
+ ... + β0 , с заданными костический полином Fçàì = p + βn-1 p
эффициентами βi .
Процедура расчета коэффициентов k1, ..., kn содержит следующие шаги.
Шаг 1. По заданным значениям чисел µ1, ..., µn находим требуемый характеристический полином замкнутой системы (4)
Fçàì = ( p - µ1 )( p - µ2 )  ( p - µn ).
Раскрывая это выражение, находим коэффициенты β0 ,  , βn-1.
Шаг 2. Находим тот же характеристический полином с помощью
формулы
Fзам = |pE – A + bk|.
Коэффициенты этого полинома будут зависеть от неизвестных
пока параметров ki.
Шаг 3. Приравнивая коэффициенты полиномов, полученных на
первом и втором шагах, получаем систему уравнений для определения неизвестных параметров ki. Решив ее, находим искомые параметры модального управления.
Замечание. В случае объектов высокого порядка объем вычислений можно сократить, записывая уравнения исходного объекта (1)
в канонической управляемой форме Фробениуса. Для объектов невысокого порядка удобнее применять описанный алгоритм.
1.3. Синтез модального регулятора
для стабилизации судна на заданном курсе
Рассмотрим судно, движущееся в плоскости (х, у). На рис. 1 изображена траектория движения судна, его вектор скорости V и текущий угол курса ϕ, требуется найти коэффициенты обратной связи
для системы управления курсом судна.
Из физических соображений ясно, что разомкнутая система
управления неустойчива, поскольку отклонение пера руля от продольной оси судна приводит к движению судна по круговой траектории. у
Линейная модель разомкнутой систеV
мы управления описывается диффеφ
ренциальным уравнением вида
φ
 + ϕ = du,
Tϕ
где Т – постоянная времени судна,
учитывающая его инерционность;
х
Рис. 1
15
d – статический коэффициент усиления канала управления, ϕ –
курсовой угол (угол между вектором скорости судна и осью х).
Переходя к изображениям по Лапласу, получаем описание канала управления в операторной форме:
Tð2 ϕ( p) + pϕ( p) = du( p).
Отсюда находим передаточную функцию
Q( p) =
d1
d
ϕ(p)
d
=
=
× 2,
u(p) Tð2 + p Tð + 1 p
где d1d2 = d.
Последней записи отвечает структурная реализация в виде последовательного соединения апериодического и интегрирующего
звеньев, показанная на рис. 2.
Здесь x2 = ϕ – угол курса; x1 = ϕ / d2 – переменная, пропорциональная угловой скорости судна; u – управляющий сигнал (угол поворота руля).
Структурная схема описывается следующей системой уравнений:
d
1
x1 (t) = - x1 (t) + 1 u(t), x2 (t) = d2 x1 (t). T
T
(5)
Корни характеристического уравнения разомкнутой системы
1
λ1 = 0, λ2 = - . Один из них равен нулю, следовательно система
Ò
асимптотически неустойчива. Для стабилизации судна на заданном
курсе введем линейную обратную связь по углу курса и угловой скорости
u(t) = -k1x1 (t) - k2 x2 (t) + k0v(t),
где v – задаваемое значение курса (управляющее воздействие).
Структура замкнутой системы управления будет иметь вид, показанный на рис. 3.
Она описывается дифференциальными уравнениями
d
d
1
x1 (t) = - x1 (t) - 1 (k1x1 (t) + k2 x2 (t)) + 1 k0v(t),
T
T
T
x2 (t) = d2 x1 (t).
(6)
В соответствии с общей методикой нужно задаться желательным
расположением корней характеристического полинома замкнутой
системы μ1 и μ2, которые должны лежать в левой полуплоскости,
16
d1
Tp +1
u
d2
x1
x2
ϕ
p
Рис. 2
v
k0
u
d1
Tp +1
x1
d2
x2
ϕ
p
– k1
–k2
Рис. 3
и рассчитать коэффициенты обратной связи k1 и k2, обеспечивающие это расположение.
Коэффициент k0 выбираем из условия равенства сигналов j и v
(угол курса равен заданному) в установившемся режиме. При этом
x1 = x2 = 0 и из уравнения (6) получаем k0 = k2.
1.4. Компьютерное моделирование
Работа выполняется в пакете MATLAB и системе SIMULINK.
При выполнении работы в командном окне MATLAB и в командном
режиме используются следующие операторы пакета MATLAB:
R=ctrb(A, B) – для вычисления матрицы управляемости (или
R=ctrb(sys));
det(R) – для вычисления определителя матрицы R;
eig(A) – для определения собственных чисел матрицы А;
y=step(A, B, C, D, i, t) – для получения графика выходного
сигнала (или y = step(sys));
k=place(A,B,M) – для получения коэффициентов обратной связи
по заданному вектору М желаемых собственных чисел замкнутой
системы;
а также команда acker и другие команды тулбокса Control пакета
MATLAB.
17
Заметим, что система в MATLAB задается в виде структуры с помощью команд ss и tf. В частности, если известны матрицы A, B, C,
D системы, то для ввода ее в MATLAB и последующего получения
графика переходной характеристики надо набрать
sys=ss(A, B, C, D),
[y, t]=step(sys); plot(t, y).
Если нужно получить не только выходной сигнал y, но и вектор
состояния Х, следует набрать [y, t, X]=step(sys).
Если не указывать выходные аргументы, а просто набрать
step(sys), то будет построен график переходной характеристики,
но никаких переменных сформировано не будет.
При моделировании разомкнутой и замкнутой систем в SIMULINK потребуются следующие блоки:
– блок передаточной функции tf;
– усилители gain;
– сумматор sum;
– генератор ступенчатого входного сигнала step;
– осциллографы Scope и XY-graph.
Они находятся в библиотеках Linear, Math, Sources, Sinks. Установка параметров этих блоков производится после двойного щелчка
мышью.
Для получения траектории судна нужно рассчитать его координаты на плоскости (х, у). Это делается по формулам
Vx = V cos ϕ (t), Vy = V sin ϕ (t),
t
x = ò Vx dt,
0
t
y = ò Vy dt.
0
Для их реализации в SIMULINK дополнительно потребуются
блоки тригонометрических функций sin, cos и интеграторы, а если
скорость судна переменная – еще и блоки произведения.
2. ЗАДАНИЕ ПО РАБОТЕ И СОДЕРЖАНИЕ ОТЧЕТА
Для заданного варианта задания требуется проанализировать
свойства разомкнутой системы (см. рис. 2), затем рассчитать коэффициенты обратной связи k1 и k2 замкнутой системы (см. рис. 3),
и исследовать зависимость величины перерегулирования Δϕ и времени переходного процесса Tp (рис. 4) от значений задаваемых собственных чисел μ1 и μ2.
18
Step Response
Step Response
x2 1.2
∆ϕ
1
0.8
0.6
Amplitude
Amplitude
0.4
0.2
0
0
5
10
Time (sec.)
Time (sec.)
t
Tpp
15
Рис. 4
При исследовании разомкнутой системы (5) необходимо:
1) проверить управляемость системы, найдя матрицу управляемости R и вычислив ее определитель;
2) выполнить анализ устойчивости, найдя собственные числа системы λ1 и λ2;
3) получить выражение для выходной переменной разомкнутой
системы и построить график ее переходной характеристики;
4) рассчитать коэффициенты k1, k2 обратной связи u(t) = – k1x1(t) –
– k2x2(t) для заданных значений μ1 и μ2;
5) построить графики переходной характеристики замкнутой системы для всех четырех случаев задания μi, указанных в вариантах;
6) выполнить моделирование в пакете MATLAB. Получить графики траектории судна на плоскости х, у при u = 1(t) в системе
SIMULINK.
Основные результаты отчета:
1. Структурная схема исходной системы и ее математическое
описание вида (5).
2. Анализ управляемости и устойчивости.
3. Формула переходной характеристики разомкнутой системы.
График этой характеристики и соответствующей траектории судна.
4. Расчет k1 и k2 и описание замкнутой системы в пространстве
состояний.
5. Графики переходной характеристики и траекторий судна для
различных значений μi.
19
6. Таблица значений ϕ, Тр и k для четырех вариантов заданий μi,
выводы об изменении Δϕ, Тр и k при изменении μi.
7. MATLAB-программы и SIMULINK-схемы для выполнения работы.
3. ПОРЯДОК ВЫПОЛНЕНИЯ
Работа выполняется в пакете MATLAB в диалоговом режиме,
а также в системе SIMULINK.
1. Получить график x2(t) и график траектории судна для разомкнутой системы.
2. Рассчитать коэффициенты k1 и k2 для заданных вариантов
значений μi.
3. Построить графики x2(t) и траектории судна замкнутой системы для различных μi при v(t) = 1(t).
4. Построить таблицу значений Δhp, Тр и ki для заданных вариантов µ i .
4. КОНТРОЛЬНЫЕ ВОПРОСЫ
1. Описать процедуру расчета коэффициентов обратной связи
при модальном управлении.
2. Получить выражение для передаточной функции замкнутой
системы (рис. 3).
3. Упростится ли вычисление ki, если:
– совпадают некоторые из собственных чисел разомкнутой
и замкнутой систем?
– совпадают некоторые из коэффициентов характеристических
полиномов разомкнутой и замкнутой систем?
4. Проверить, изменяется ли числитель передаточной функции
системы при модальном управлении.
5. Как изменятся свойства системы, если одно из собственных
чисел замкнутой системы окажется совпадающим с нулем ее передаточной функции?
6. Найти коэффициенты ki и передаточную функцию замкнутой
системы, если µ1 = µ2 = 0.
20
5. ВАРИАНТЫ ЗАДАНИЙ
№
1
2
3
4
5
6
7
8
9
10
d1
1
1
1
1,5
1,5
1,5
1
1
1
2
d2
1
1
1
1
1
1
1,5
1,5
1,5
1
T
1
1,5
2
2,5
3
3,5
4
4,5
5
5,5
m
2
2,5
3
3,5
4
4,5
5
5,5
6
5,5
№
11
12
13
14
15
16
17
18
19
20
d1
2
2
1
1
1
2,5
2,5
2,5
4
3
d2
1
1
2
2
2
1,5
1,5
1,5
4
3
T
6
5,5
5
4,5
4
3,5
3
2,5
2
1,5
m
5
4,5
4
3,5
3
2,5
2
1
2
1,5
В каждом из вариантов рассмотреть четыре случая задания корней характеристического полинома замкнутой системы:
а) μ1= μ2= – m;
б) μ1= – m, μ2=10 μ1;
в) µ1,2 = α ± jβ, α = -m, β = 0,1α;
г) µ1,2 = α ± jβ, α = -m, β = α.
Расчеты производить для нулевых начальных условий.
21
Лабораторная работа № 3
ПРОГРАММНОЕ УПРАВЛЕНИЕ
ПОДВИЖНЫМИ ОБЪЕКТАМИ
Цель работы: освоить методику расчета и компьютерного моделирования программного управления, обеспечивающего перевод
подвижного объекта в заданное конечное положение при различных ограничениях на управляющие воздействия.
1. ТЕОРЕТИЧЕСКИЕ СВЕДЕНИЯ
1.1. Принцип программного управления
Одна из классических задач теории управления – это задача перевода подвижного объекта (динамической системы) из заданного
исходного состояния X(0) в заданное конечное состояние X(T) за
фиксированное время Т. Она известна как задача терминального
управления с фиксированным временем.
Существуют два подхода к ее решению. Первый из них использует программное управление, второй – управление с помощью обратной связи. В данной работе рассматривается первый подход, согласно которому управляющее воздействие u(t) рассчитывается заранее
и формируется с помощью специального генератора (рис. 1).
Будем считать, что управляемая динамическая система описана
матричным уравнением в пространстве состояний:
 = AX + bu, X(0) = X , X
(1)
0
где X ∈
столбец; u – скалярный управляющий сигнал.
Цель управления – перевести систему в заданное состояние
X(T). Чтобы получить уравнения для определения сигнала u(t),
обеспечивающего достижение этой цели, воспользуемся интеRn – вектор состояния; А – квадратная матрица; b – вектор-
ГЕНЕРАТОР
u
CИСТЕМА
Рис. 1
22
X
гралом свертки, связывающим начальное и конечное состояния
системы
T
X(T) = e ÀT X(0) + ò Q(T - t)u(t)dt, (2)
0
где Q(t) = eAt b – импульсная весовая функция системы; eAt – матричная экспонента.
Функция Q(t) представляет собой вектор-столбец, элементами
которого являются известные скалярные функции qi (t).
Введем обозначения M = X(T) – eATX(0), F(t) = Q(T – t) и воспользуемся стандартной формулой для скалярного произведения двух
функций
T
(f, g) = ò f (t) g(t)dt.
0
Это позволяет переписать матричное равенство (2) в виде совокупности n скалярных равенств
(f1 , u) = m1,
(3)
(fn , u) = mn.
Здесь mi – компоненты постоянного вектора М; fi – компоненты
вектор-функции F(t).
С алгебраической точки зрения равенства (3) представляют собой систему уравнений для определения неизвестной функции u(t).
T
В математике интегралы вида
ò fi (t)u(t)dt
называют моментами
0
функции u(t), поэтому полученная задача состоит в восстановлении
функции по ее моментам.
С геометрической точки зрения функции fi и u удобно рассматривать как векторы (они принадлежат гильбертову пространству).
Тогда числа mi представляют собой проекции вектора u на направления, задаваемые векторами f1, ..., fn. Таким образом, геометрически задача сводится к определению вектора по его проекциям. Поскольку в нашем случае число проекций меньше размерности вектора, задача имеет бесчисленное множество решений (существует
много векторов, у которых часть компонент совпадает, а остальные
различаются).
На практике обычно бывает нужно найти конкретное управление, удовлетворяющее уравнениям (3), т.е. выбрать из множества
23
возможных решений какое-то одно. Для этого осуществим параметризацию решения, а именно будем искать его в виде известной
функции, зависящей от n параметров, которые необходимо определить. Рассмотрим несколько простых вариантов параметризации.
1.2. Кусочно-постоянное управление
Кусочно-постоянное управление с фиксированными моментами
переключения. Разобьем интервал управления на n равных участков и обозначим амплитуду управляющего сигнала на каждом из
участков через c1, ..., cn (рис. 2).
Тем самым от произвольных функций u(t) мы перешли к более
узкому классу кусочно-постоянных функций, характеризуемых n
параметрами c1, ..., cn.
Для таких функций уравнения (3) преобразуются к виду
t1
T
0
tn -1
c1 ò fi dt + ... + cn ò fi dt = mi ;
i = 1, ..., n.
В результате получается система n линейных алгебраических
уравнений для определения n неизвестных параметров c1, ..., cn.
Кусочно-постоянное управление с постоянной амплитудой
(релейное управление). Потребуем, чтобы управление принимало
только два значения ±с (так называемое релейное или сигнатурное
управление), а в качестве неизвестных параметров возьмем моменты переключения t1, t2, ..., tn–1. Примерный вид такого управления
показан на рис. 3.
В этом случае система уравнений (3) принимает вид
t1
c ò f1dt ± ... ± c
0
T
ò
f1dt = mi ; i = 1,..., n.
tn-1
После вычисления интегралов получаем систему из n нелинейных алгебраических уравнений с n неизвестными t1, ..., tn-1, c.
с1
u
u
...
с2
сn
...
с
t1
T
t
0
tn–1
T
t1
t2
tn–1
–с
Рис. 2
24
Рис. 3
t
Решить такую систему значительно сложнее, чем в предыдущем
случае, однако получаемое управление обладает двумя замечательными свойствами.
Во-первых, это управление будет иметь минимально возможную
амплитуду среди всех возможных управлений. Другими словами
полученное управление будет обладать минимальной чебышевской
нормой
u(t) ¥ = max u(t) ® min.
0£t£T
Во-вторых, для систем n-го порядка с вещественными корнями
характеристического полинома число переключений будет не более, чем n – 1. Это утверждение известно как теорема Фельдбаума
об n интервалах.
1.3. Управление с минимальной энергией
Недостатком двух предыдущих вариантов параметризации было наличие скачков в управляющем сигнале, а также повышенные
энергетические затраты на управление. Чтобы избежать этих недостатков, поставим задачу отыскания управления с минимальной
энергией.
Энергию управляющего сигнала будем оценивать по формуле
T
E = u(t) 22 = (u(t), u(t)) = ò u2dt. (4)
0
Рассматривая величину Е как минимизируемый функционал,
а равенства (3) как ограничения, получаем задачу на условный экстремум.
Для ее решения строим функцию Лангранжа
L = (u, u) + l1((f1, u) – m1) + ... + ln((f1, u) – mn),
дифференцируем ее по u и приравниваем производную нулю:
Lu¢ = 2u + λ1f1 + ... + λn fn = 0.
Отсюда получаем следующее выражение для функции u:
u(t) = c1 f1(t) + ... + cn fn(t),
(5)
где сi = –li/2 – постоянные коэффициенты, подлежащие определению.
Таким образом, оптимальное управление, обладающее минимальной энергией, должно быть линейной комбинацией весовых
25
функций q1(t), ..., qn(t) управляемой системы (1), взятых в обратном
времени.
Подставляя выражение (5) в формулы (3), получаем систему линейных алгебраических уравнений относительно неизвестных коэффициентов c1, ..., cn:
é(f1 , f1 ) ... (f1 , fn ) ù é c1 ù ém1 ù
ê
úê ú ê ú
ê.....................ú ê...ú = ê... ú .
ê
úê ú ê ú
ê
úê ú ê ú
ë(fn , f1 ) ... (fn , fn ) û ë c û ëmn û
Ее удобно записать в матричной форме: WC = M, где W =
n
T
W = ò F (t) F T (t)dt.
0
Отсюда находим окончательное выражение для искомого управления:
(6)
u(t) = FT(t)C = FT(t)W–1M.
В более подробной записи формула для матрицы W примет вид
T
T
0
0
T
W = ò F (t)F T (t)dt = ò e At bbT e A t dt.
(7)
Симметричная матрица W, определяемая этой формулой, называется грамианом управляемости системы (1). Она играет важную
роль в теории управления.
Полученное управление (6) отличается от всех других тем, что
оно имеет минимально возможную энергию (4).
Найдем величину этой энергии:
T
T
0
0
E = ò (F T W-1M)T (F T W-1M)dt = MT W-1 (ò FF T dt)W-1M = MT W-1M. (8)
Например, если положить X(0) = 0; X(T) = Xк (задача о переводе системы из начала координат в заданное состояние за время Т),
то энергия входного сигнала определяется квадратичной формой
Е = ХкТW–1Xк.
1.4. Задача об управлении движением двух материальных тел
Рассмотрим задачу о встрече двух материальных тел А1 и А2,
движущихся вдоль прямой (рис. 4). Например, речь может идти
о движении двух самолетов, о движении двух автомобилей и т. д.
Дополнительное условие состоит в том, что в момент встречи оба
26
А1
0
u
v
y0/2
А2
y0
x, y
Рис. 4
тела должны иметь одинаковую скорость (это необходимо, например, при заправке самолета в воздухе). Величина этой скорости, а
также место и время встречи задаются заранее.
Уравнения движения тел получаются из второго закона Ньютона F = ma. В простейшем случае движения двух тел с единичной
массой при отсутствии сопротивления среды получаем дифференциальные уравнения второго порядка
 = u, x(0) = x0 , x (0) = x 0 ;
x

y = v, y(0) = y0 , y (0) = y0 .
 , y , x , y , x, y – ускорения, скорости и положения тел;
Здесь x
u и v – управляющие силы, например, силы тяги двигателей.
Будем считать, что в начальный момент времени тела были неподвижны и имели координаты 0 и y0 соответственно. Потребуем,
чтобы через заданное время Т оба тела имели одно и то же положение x(T) = y(T) = y0/2 и одинаковые скорости x (T) = y (T) = 1. Нужно
рассчитать три варианта управлений u(t) и v(t), обеспечивающих достижение этой цели.
Примерный вид траекторий на фазовой плоскости, характеризующих желаемое движение тел, показан на рис. 5. На нем точки А1
и А2 соответствуют начальному состоянию управляемых тел, точка А – конечному состоянию.
Ясно, что вид траекторий будет зависеть от типа закона управления. Для отыскания требуемых управлений удобно перейти к матричной форме описания системы (1). В данном случае оба объекта
описываются одинаковыми матрицами A, b вида
é 0 1ù
ú,
A=ê
êë0 0úû
с граничными условиями
é 0ù
b = ê ú, êë1úû
(9)
é y0 /2ù
é 0ù
é y0 ù
ú.
X(0) = ê ú , Y (0) = ê ú ; X(T) = Y (T) = ê
êë 1 úû
êë0úû
êë 0 úû
27
. .
x, y
A
1
A1
A2
y0/2
0
x, y
y0
Рис. 5
Отсюда для вектор-функций Q(t) и F(t) получаем (покажите это!)
étù
éT - t ù
ú.
Q(t) = ê ú ,
F (t) = ê
êë1úû
êë 1 úû
Подставляя эти формулы в выражения (3) и выполняя три рассмотренные выше варианта параметризации, можно получить разные варианты управлений для рассматриваемых объектов.
В качестве примера найдем управления, обеспечивающие перевод тела A1 из начала координат в состояние A с координатами (1, 1)
за время Т = 1.
Уравнения (3) в этом случае имеют вид (1 – t, u) = 1; (1, u) = 1. Заменяя для удобства первое из этих двух уравнений их разностью,
получаем (t, u) = 0; (1, u) = 1 или в интегральной форме
1
ò tudt = 0;
0
1
ò udt = 1. (10)
0
При первом варианте параметризации (кусочно-постоянное
управление) функция u = c1 на интервале 0 ≤ t < 1/2 и u = c2 на интервале 1/2 < t ≤ 1. В этом случае из уравнений (10) находим
c1
t2
2
0,5
0
+ c2
t2
2
1
0,5
= 0;
0,5c1 + 0,5c2 = 1,
откуда с1 = 3, с2 = –1, т. е. управление имеет вид, показанный на
рис. 6, а.
28
При втором варианте параметризации (релейное управление)
функция u = c на интервале 0 ≤ t < t1 и u = –c на интервале t1 < t < 1,
где t1 – неизвестный момент времени. В данном случае из уравнений (10) получаем
ct12/2 – c(1 – t12)/2 = 0, ct1 – c(1 – t1) = 1,
откуда t1 = 2 / 2 » 0,707; c = 1 + 2 » 2,41.
График соответствующего управления приведен на рис. 6, б.
При третьем варианте параметризации (управление с минимальной энергией) управление имеет вид u = c1 + c2t. Из системы (10) получаем
1
1
t2
t3
t2
(c1 + c2 ) = 0; (c1t + c2 ) = 1,
2
3 0
2 0
или c1/2 + c2/3 = 0; c1 + c2/2 =1. Отсюда c1 = 4; c2 = –6, u = 4 – 6t.
График этого управления пересекает ось абсцисс в точке t = 2/3
(рис. 6, в). Зависимость скорости и положения тела от времени характеризуются формулами x2 = 4t – 3t2; x1 = 2t2 – t3.
Вид траекторий на фазовой плоскости (x, x ) для трех рассмотренных примеров приведен на рис. 7, а–в. Все траектории выходят
из начала координат и заканчиваются в точке с координатами (1, 1).
а)
в)
б)
u
4
u
2
1
0
–1
0,5
1 t
u
2
2,41
0,5
0,707 t
1
–2
–2,41
0,667 t
0,5
1
Рис. 6
а)
а)
б)
1.6
в)
1.8
1.6
1.2
0.8
1.2
0.8
0.8
0.4
0
0.4
0.4
0
0.4
0.8 1
0
1.4
1.2
0
0.4
0.8 1
0
0
0.4
0.8
1
Рис. 7
29
2. ЗАДАНИЕ ПО РАБОТЕ
И СОДЕРЖАНИЕ ОТЧЕТА
1. Составить систему уравнений (3) для материальных тел А1
и А2 (см. рис. 4) в соответствии с заданными значениями y0 и Т. Найти ее решение для каждого из трех типов параметризации управляющего сигнала.
2. Для каждого из найденных управлений построить графики
движения тел в зависимости от времени.
3. Нарисовать траектории движения в фазовой плоскости.
4. Привести схемы моделирования для выполнения работы в пакете SIMULINK. Схема для каждого тела должна содержать два интегратора, генератор управляющего сигнала и осциллографы Scope
и XY-graph для регистрации результатов. Генератор кусочно-постоянного управляющего сигнала может быть построен с помощью блоков
step и сумматора, либо с помощью блока Signal Builder из библиотеки Sources. Для получения линейно-изменяющегося управляющего сигнала следует использовать блок Ramp из той же библиотеки.
Он имеет параметры Slope и Initial output для задания наклона и начальной точки прямой.
5. Привести программы для выполнения работы в MATLAB. Они
должны включать формирование управляющих сигналов по формулам, полученным в отчете, создание объектов sys1 и sys2 с помощью
конструктора ss и расчет выходных сигналов объектов с помощью
команды lsim, а также построение всех необходимых графиков.
Для получения кусочно-постоянных сигналов удобно использовать
знаки неравенств и логические операторы. Например, команды
t=1:10; u=t>3&t<7
сформируют последовательность
u = 0 0 0 1 1 1 0 0 0 0
3. ПОРЯДОК ВЫПОЛНЕНИЯ
Работа выполняется в MATLAB и SIMULINK.
1. При выполнении следует получить графики координат, скоростей и фазовые траектории для всех вариантов управлений. При
этом управляющие сигналы, полученные в отчете, подаются на модели объектов (sys1 и sys2 в MATLAB или схемы в SIMULINK) и регистрируются выходные сигналы.
30
2. Для оценки точности решения задачи о встрече двух тел сформировать функционал J(t) = (y - x)2 + (y - x )2 и проанализировать
его поведение во времени для разных законов управления.
3. Определить нормы управляющих воздействий u 1,2,¥ для
каждого из вариантов управлений. Полученные данные оформить
в виде таблицы норм.
4. Исследовать чувствительность системы при разных законах
управления к трем типам внешних помех:
a) наличие постоянного ветра (u ± 0,1, v ± 0,1);
 + 0,1x = u, y + 0,1y = v);
б) наличие сопротивления воздуха (x
в) уменьшение мощности двигателей до величин 0,9u, 0,9v.
Полученные данные оформить в виде таблицы для функционала J(T).
4. КОНТРОЛЬНЫЕ ВОПРОСЫ
1. Найти управления для перевода объекта из точки А в начало координат (см. рис. 4) и построить графики управлений с минимальной энергией и с минимальной амплитудой.
2. Найти управление для перевода объекта из точки А в точку А2
и построить графики управлений с минимальной энергией и с минимальной амплитудой.
3. В начальный момент оба объекта находятся в одной точке пространства и имеют скорости 1 и 2. Найти управление для перевода
их за одну секунду в эту же пространственную точку с нулевыми
конечными скоростями.
4. Используя формулу (7), найти грамиан управляемости для исследуемого объекта (9). Рассчитать управление для своего варианта
по формуле (6) и сравнить его с полученным в отчете.
5. Рассчитать энергию управления для каждого из объектов по
формуле (9) и сравнить ее с полученной в отчете.
6. Как изменятся решения, полученные в отчете, если принять
массу второго тела равной 0,5? Нарисовать траектории в фазовой
плоскости.
7. Оценить чувствительность системы при разных законах
управления к трем типам внешних помех, описанных в п. 3.4.
8. Рассмотреть аналогичную задачу о встрече двух тел на плоскости.
31
5. ВАРИАНТЫ ЗАДАНИЙ
Во всех вариантах принять
x(0) = 0, y(0) = y0 , x (0) = y (0) = 0; x(T) = y(T) = y0/2, x (T) = y (T) = 1.
32
№
1
2
3
4
5
6
7
8
9
10
y0
1
2
10
8
3
2
1
8
2
18
T
2
1,5
5
4
1
3
0,7
2
5
15
№
11
12
13
14
15
16
17
18
19
20
y0
14
18
4
6
2
3
7
2
5
15
T
7
9
1
3
6
2
2
4
2
4
Лабораторная работа № 4
ЭКСТРЕМАЛЬНЫЕ НАЧАЛЬНЫЕ УСЛОВИЯ
ЛИНЕЙНОЙ СИСТЕМЫ
Цель работы: освоить методику расчета начальных условий линейной динамической системы для получения выходного сигнала
с экстремальной энергией, используя пакет MATLAB.
1. ТЕОРЕТИЧЕСКИЕ СВЕДЕНИЯ
1.1. Задача об экстремуме квадратичной формы на сфере
Квадратичной формой от переменных x1, , xn называется выражение вида
n
Q(X) = å wij xi xj . (1)
i, j
Одна из классических задач оптимизации состоит в том, чтобы
найти экстремумы этой квадратичной формы при условии
x12 +  + xn2 = 1. (2)
Геометрически этому уравнению отвечает единичная сфера в n-мерном пространстве.
Матричная запись уравнений (1), (2) имеет вид
Q(X) = XT WX, XT X = 1,
где Х – вектор с компонентами xi; W – симметричная квадратная
матрица с элементами wij; T – символ транспонирования.
Для отыскания экстремума формируем составной критерий
L = XT WX + λ(XT X -1),
где λ – множитель Лагранжа.
Приравнивая производную по X нулю, получаем
2WX - 2λX = 0,
WX = λX.
Последнее равенство означает, что искомые экстремальные направления совпадают с собственными векторами H1, ..., Hn матрицы W, которые вследствие симметричности матрицы W образуют
ортонормальную систему.
33
Вычисляя значения квадратичной формы (1) при X = Hi, получаем Q(H1) = λ1, ..., Q(Hn) = λn, т. е. собственные числа матрицы W
численно равны экстремальным значениям квадратичной формы
на единичной сфере. Максимальное и минимальное из них задают
границы изменения квадратичной формы при X = 1:
λn £ Q(X) £ λ1.
Задача о максимуме квадратичной формы на сфере имеет многочисленные приложения в технических дисциплинах.
1.2. Задача об экстремальных начальных условиях
Пусть исследуемая линейная динамическая система задана матричным описанием вида
 = AX, y = cX,
X
где Х – вектор с компонентами x1, , xn ; А – квадратная матрица;
с – вектор-строка с элементами c1, , cn .
Требуется найти вектор начальных условий Х0 единичной длины X0 = 1, при котором энергия сигнала y(t), определяемая фор¥
мулой Ey = ò y2 (t)dt, была бы максимальной.
0
Рассмотрим два пути решения этой задачи.
Как известно, выходной сигнал системы при произвольных начальных условий может быть записан в виде
y = c1e p1t +  + cn e pnt ,
где pi – собственные числа матрицы А; ci – константы, зависящие от
вектора начальных условий Х0.
Энергия сигнала у определяется формулой
¥
¥
Ey = ò y2dt = ò (c1e p1t + + cn e pnt )2 dt. (3)
0
0
Первый путь, приводящий к громоздким вычислениям, состоит
в том, чтобы явно выразить константы ci через начальные условия,
подставить их в формулу (3) далее искать экстремум полученного
выражения.
Второй путь основан на выражении сигналов X(t) и y(t) через матричную экспоненту
34
X(t) = e At X0 ,
y(t) = ce At X0 . (4)
Формально матричная экспонента определяется с помощью бесконечного ряда
e At = E + At +
( At)2 ( At)3
+
+ ,
2!
3!
аналогичному ряду Тейлора для обычной экспоненты.
На практике столбцы матричной экспоненты можно получить
 = AX при начальных успутем решения системы уравнений X
T
ловиях, равных столбцам единичной матрицы éë1 0  0ùû , , éë0 0 
T
T
¥
 0ùû , , éë0 0  1ùû .
Перепишем формулу (3) для энергии в виде Ey = ò yT ydt и под0
ставим вместо у(t) его аналитическое выражение (4). Выполняя
несложные преобразования, получим
æ¥
ö
AT t T At ÷÷÷
T
e
c
ce
d
t
÷÷X0 = X0 W0 X0 . ò
÷ø
çè 0
ç
Ey = X0T çç
ç
(5)
Обозначим через W0 выражение, заключенное в скобки.
¥
T
W0 = ò e A t cT ce At dt. (6)
0
Оно представляет собой симметричную числовую матрицу, которая в теории управления называется грамианом наблюдаемости.
Тем самым задача свелась к исследованию квадратичной формы X0T W0 X0 при ограничении X0T X0 = 1, т. е. задача об экстремуме квадратичной формы на сфере? решение которой описано выше
(п. 1.1).
1.3. Пример расчета экстремальных начальных условий
Рассмотрим задачу об оптимальных начальных условиях для системы второго порядка, описываемой дифференциальными уравнениями
x = -x,
(7)
y = x - 0,5y.
Перейдем к матричной форме записи:
é xù
é-1
0 ù
̂
ú , с = [0 1].
X = AX; y = cX, где X = ê ú , A = ê
ê yú
ê 1 -0,5ú
ë û
ë
û
(8)
35
y
1
y0
Нужно исследовать зависимость энергии выходного сигнала
X0
x0 1
Рис. 1
¥
Ey = ò y2 (t)dt
х
0
от начальных условий x0, y0, ответив на следующие вопросы:
а) в каких пределах может изменяться энергия Eу, если векторы начальных условий имеют
единичную длину x2 + y2 = 1, т. е. стартовые точки лежат на единичной окружности в фазовой плоскости (рис. 1)?
б) при каких начальных условиях энергия Eу будет наибольшей?
в) как расположены на плоскости (x, y) точки начальных условий, для которых Eу = 1?
Для получения ответа на эти вопросы надо найти матричную
экспоненту, вычислить грамиан наблюдаемости и найти его собственные векторы.
Решая дифференциальные уравнения (7) при начальных условиях x0 = 1, y0 = 0, получаем
x = e-t , y = 2e-0,5t - 2e-t .
Здесь показатели экспонент – это собственные числа матрицы А
p1 = -1, p2 = -0,5.
Решая те же уравнения при x0 = 0, y0 = 1, получаем
x = 0,
y = e-0,5t .
Следовательно матричная экспонента имеет вид
é
e-t
e At = êê
êë-2e-t + 2e-0,5t
ù
ú.
-0,5t úú
e
û
0
Умножаем ее на вектор-строку c = [0 1] :
ce At = éëê-2e-t + 2e-0,5t e-0,5t ùûú .
Для вычисления грамиана наблюдаемости используем формулу (4)
¥é
-2(e-t + 2e-0,5t úù é
-t
-0,5t
W0 = ò êê
ú êë-2e + 2e
0
,
5
t
ê
úû
e
0 ë
36
1 é2 2ùú
e-0,5t ùú dt = ê
.
û
3 ëê2 3ûú
Другой способ нахождения грамиана наблюдаемости основан на
решении матричного уравнения Ляпунова
W0 A + A T W0 = -cT c.
В скалярной записи оно представляет собой систему четырех
уравнений с четырьмя неизвестными – элементами грамиана наблюдаемости. Его можно также использовать для проверки правильности найденной матрицы W0.
Далее надо найти собственные числа и собственные вектора грамиана наблюдаемости. Выписываем характеристический полином
грамиана и находим его корни
5
2
λE - W0 = λ2 - λ + ;
λ1 = 1,52; λ2 = 0,146.
3
9
Соответствующие собственные вектора имеют вид
é0,6154ù
é 0,7882 ù
ú,
ú.
H1 = ê
H2 = ê
ê
ú
ê
ú
ë 0,7882û
ë-0,6154û
Они удовлетворяют уравнениям
W0 H1 = λ1 H1, W0 H2 = λ2 H2 ,
H1 = H2 = 1,
и могут быть найдены в пакете MATLAB с помощью команды
[H, L]=eig(W0).
Числа λ1, λ2 характеризуют максимальную и минимальную
энергию Еу, которую можно получить при начальных условиях,
взятых с единичной окружности.
Начальные условия, при которых достигаются экстремальные
значения энергии, определяются собственными векторами H1, H2 .
Экстремальные кривые уМ, уm показаны на рис. 2 слева. Справа показаны соответствующие траектории на плоскости (x, y).
1
0.8
0.5
0.6
уМ
y
0.4
y
0.2
уm
0
-0.5
0
-0.2
0
5
10
t
15
Рис. 2
-1
-1
-0.5
0
0.5
1
x
37
Они начинаются на единичной окружности и заканчиваются в начале координат.
Графики построены в MATLAB с помощью команд
sys=ss(A,[],c,0);
t=0:.01:15;
y1=initial(sys,H(:,1),t);
y2=initial(sys,H(:,2),t);
plot(t,[y1,y2]),grid
plot (y1,y2),grid
Рассчитаем отношение максимальной и минимальной энергий
λ1 / λ2 = 10,4.
Следовательно, при одной и той же длине вектора начальных условий, но разных его направлениях, энергия выходного сигнала может изменяться более, чем в 10 раз.
Мы получили ответ на все поставленные вопросы. Диапазон изменения энергии Ey при условии x02 + y02 = 1, определяется двусторонним неравенством λ2 £ Ey £ λ1, где λ2 и λ2 – максимальное и
минимальное собственные числа грамиана наблюдаемости.
Максимально возможная энергия колебаний равна λ1. Она достигается, если вектор начальных условий совпадает по направлению с первым собственным вектором грамиана: X0 = H1 . Минимальная энергия получается, если в качестве начальных условий
взять второй собственный вектор X0 = H2 .
Множество начальных условий, для которых энергия колебаний
будет единичной, задается уравнением X0T W0 X0 = 1. В теории
управления этот эллипс называется эллипсом наблюдаемости. Чем
сильнее он вытянут вдоль некоторого направления, тем хуже оно
наблюдаемо (для получения единичной энергии на выходе надо задавать большие начальные условия).
2. ЗАДАНИЕ ПО РАБОТЕ И СОДЕРЖАНИЕ ОТЧЕТА
В лабораторной работе исследуется система, представляющей
собой последовательное соединение двух апериодических звеньев
(рис. 3).
Для заданных значений а1, а2 требуется:
1) привести описание системы к матричной форме вида (8):
é xù
 = AX; y = cX, X(0) = X ; X = ê ú . X
(9)
0
ê yú
ë û
Найти собственные числа матрицы А. Заметим, что в нашем случае c = éë0 1ùû , а матрица А треугольная;
38
x0
1
p + a1
y0
x
1
p + a2
y
Рис. 3
2) решить систему (9) для начальных условий вида [1 0]T и [0 1]T
и построить графики х(t), y(t);
3) найти матричную экспоненту и вычислить грамиан наблюдаемости системы W0;
4) найти собственные числа и собственные векторы Н1, Н2 грамиана W0 и оценить диапазон возможного изменения энергии сигнала у(t);
5) составить MATLAB-программу для численного решения задачи, путем просмотра с шагом 3° точек единичной окружности, выбирая их в качестве начальных условий.
3. ПОРЯДОК ВЫПОЛНЕНИЯ
1. Ввести матрицы А, c, и с сформировать помощью команды ss
систему sys. С помощью команды eig найти собственные числа матрицы А. Используя команду initial для начальных условий [1; 0]
и [0; 1], получить графики сигналов Х1(t), Х2(t),
2. Используя команды gram, eig, найти грамиан W0 и его собственные векторы Н1, Н2.
3. Получить графики выходных сигналов для Х0 = Н1, Х0 = Н2.
С помощью команды trapz найти энергию этих сигналов и сравнить
ее с собственными числами грамиана.
4. Используя оператор цикла for, выполнить п. 5 задания. Сравнить полученные экстремальные начальные условия с собственными векторами Н1, Н2.
5. Получить фазовые траектории для всех случаев моделирования исходной и преобразованной систем. Используя функцию
ezplot, построить эллипс наблюдаемости.
39
4. КОНТРОЛЬНЫЕ ВОПРОСЫ
x0
1. Записать квадратичную форму с матрицей А из своего варианта и найти ее экстремумы на единичной окружности.
2. Найти матрицы А и с для схемы рис. 4.
Найти матричную экспоненту, взяв a1, a2 из
своего варианта.
3. Найти грамиан наблюдаемости системы из своего варианта, заменив матрицу
с на следующую:
x1
1
p + a1
y
x2
1
p + a2
Рис. 4
a) c = [1 1],
á) c = [ 1
0 ],
â) c = [1 2],
ã) c= [ 1 – 1].
5. ВАРИАНТЫ ЗАДАНИЙ
Исследуемая система представляет собой последовательное соединение двух апериодических звеньев (рис. 3).
Значения параметров а1 и а2 для каждого из вариантов приведены в таблице.
№ 1
2
3
4
5
6
7
8
9 10 11 12 13 14 15 16 17 18 19 20
а1 1
1
1
1
2
2
2
3
3
4
2
3
4
5
3
4
5
4
5
5
а2 2
3
4
5
3
4
5
4
5
5
1
1
1
1
2
2
2
3
3
4
40
Содержание
Лабораторная работа № 1. Синтез пид-регулятора........... 1. Теоретические сведения......................................... 2. Задание по работе и содержание отчета.................... 3. Порядок выполнения............................................ 4. Контрольные вопросы........................................... 5. Варианты заданий................................................ Лабораторная работа № 2. Модальное управление............. 1. Теоретические сведения......................................... 2. Задание по работе и содержание отчета.................... 3. Порядок выполнения ............................................ 4. Контрольные вопросы........................................... 5. Варианты заданий................................................ Лабораторная работа № 3. Программное управление подвижными объектами.................................................... 1. Теоретические сведения......................................... 2. Задание по работе и содержание отчета.................... 3. Порядок выполнения............................................ 4. Контрольные вопросы........................................... 5. Варианты заданий................................................ Лабораторная работа № 4. Экстремальные начальные
условия линейной системы............................................ 1. Теоретические сведения......................................... 2. Задание по работе и содержание отчета.................... 3. Порядок выполнения............................................ 4. Контрольные вопросы........................................... 5. Варианты заданий................................................ 3
3
11
11
11
12
13
13
18
20
20
21
22
22
30
30
31
32
33
33
38
39
40
40
41
Для заметок
42
Для заметок
43
Документ
Категория
Без категории
Просмотров
0
Размер файла
2 095 Кб
Теги
mironovskii
1/--страниц
Пожаловаться на содержимое документа