МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ
РОССИЙСКОЙ ФЕДЕРАЦИИ
ФЕДЕРАЛЬНОЕ АГЕНТСТВО ПО ОБРАЗОВАНИЮ
РОССИЙСКОЙ ФЕДЕРАЦИИ
Государственное образовательное учреждение
высшего профессионального образования
САНКТ-ПЕТЕРБУРГСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ
АЭРОКОСМИЧЕСКОГО ПРИБОРОСТРОЕНИЯ
А.И. Панферов, А.В. Лопарев
КОМПЬЮТЕРНЫЙ АНАЛИЗ И СИНТЕЗ СИСТЕМ
ОРИЕНТАЦИИ, СТАБИЛИЗАЦИИ И НАВИГАЦИИ
Учебное пособие
Санкт-Петербург
2008
Учебное пособие содержит описание лабораторных работ по исследованию динамических и точностных свойств систем ориентации, стабилизации и
навигации. Каждая работа содержит методические указания к выполнению заданий, описание программы моделирования, контрольные вопросы. Приводятся необходимые сведения из теории фильтрации и теории управления. В качестве среды моделирования используется пакет Matlab 6.5.
Пособие предназначено для студентов технических вузов, обучающихся
по специальности 160402 и направлению 200100.
Рецензенты:
кафедра информационно-навигационных систем Санкт-Петербургского
государственного университета информационных технологий, механики и
оптики;
кандидат технических наук А.Е. Елисеенков (ГНЦ РФ – ЦНИИ «Электроприбор»)
Тираж 100 экз.
2
ОГЛАВЛЕНИЕ
ПРЕДИСЛОВИЕ.......................................................................................................... 5
1.
СИНТЕЗ СИСТЕМЫ СТАБИЛИЗАЦИИ СУДНА ПО КУРСУ С
ЗАДАННЫМ КАЧЕСТВОМ ПЕРЕХОДНОГО ПРОЦЕССА ................................ 7
1.1. Методические указания по подготовке к работе..................................... 7
1.2. Описание программы моделирования .................................................... 14
1.3. Порядок выполнения лабораторной работы .......................................... 16
1.4. Оформление отчета ................................................................................... 22
1.5. Контрольные вопросы .............................................................................. 22
2.
ИССЛЕДОВАНИЕ ТОЧНОСТИ ОПРЕДЕЛЕНИЯ УГЛОВ ОРИЕНТАЦИИ
ПОДВИЖНОГО ОБЪЕКТА С ИСПОЛЬЗОВАНИЕМ БЕСПЛАТФОРМЕННОЙ
ИНЕРЦИАЛЬНОЙ ВЕРТИКАЛИ ........................................................................... 24
2.1. Методические указания по подготовке к работе................................... 24
2.2. Описание программы моделирования .................................................... 33
2.3. Порядок выполнения лабораторной работы .......................................... 36
2.4. Оформление отчета ................................................................................... 43
2.5. Контрольные вопросы .............................................................................. 43
3.
СИНТЕЗ РЕГУЛЯТОРА ДЛЯ СИСТЕМЫ УПРАВЛЕНИЯ УПРУГИМ
ОБЪЕКТОМ ............................................................................................................... 45
3.1. Методические указания по подготовке к работе................................... 45
3.2. Описание программы моделирования .................................................... 53
3.3. Порядок выполнения лабораторной работы .......................................... 54
3.4. Оформление отчета ................................................................................... 59
3.5. Контрольные вопросы .............................................................................. 59
4.
СИНТЕЗ АЛГОРИТМА ЭКСТРАПОЛЯЦИИ УЗКОПОЛОСНОГО
СЛУЧАЙНОГО СИГНАЛА ..................................................................................... 61
4.1. Методические указания по подготовке к работе................................... 61
4.2. Описание программы моделирования .................................................... 66
3
4.3. Порядок выполнения лабораторной работы .......................................... 68
4.4. Оформление отчета ................................................................................... 71
4.5. Контрольные вопросы .............................................................................. 71
ЛИТЕРАТУРА ........................................................................................................... 73
ПРИЛОЖЕНИЯ ........................................................................................................ 75
4
ПРЕДИСЛОВИЕ
Развитие средств вычислительной техники предоставляет все больше
возможностей при решении задач теории автоматического управления. Вычислительные устройства используются как на этапе синтеза при предварительном
расчете характеристик системы, так и непосредственно в контуре управления
при ее практической реализации.
Одним из наиболее популярных инструментов научных исследований и
математических расчетов стала система Matlab, построенная на расширенном
представлении и применении матричных операций. Матрицы широко применяются при построении и математическом моделировании динамических систем и объектов, они являются основой составления и решения уравнений состояния. Однако в настоящее время Matlab далеко вышла за пределы специализированной матричной системы и является одним из наиболее мощных универсальных средств компьютерных расчетов и моделирования. При этом большое
внимание уделяется наглядности и удобству в использовании системы, что
наглядно проявляется в пакете расширения Matlab – Simulink.
Обладая широкими возможностями о созданию различных демонстрационных примеров, система Matlab нашла свое применение и в учебном процессе.
Изданная учебно-методическая литература по анализу и синтезу динамических
систем, так или иначе опирающаяся на Matlab [1-6], охватывает широких
спектр теоретических и практических вопросов, однако исследовательскому
аспекту в подобных работах отведена второстепенная роль. Вместе с тем, подготовка специалистов и магистров по направлению 652300 «Системы управления движением и навигация» в высших учебных заведениях подразумевает
приобретение студентами практических навыков проектирования сложных динамических систем. Такие навыки могут быть получены в ходе выполнения лабораторных работ с применением средств вычислительной техники, и предла-
5
гаемое учебное пособие может восполнить недостаток учебного материала в
данной области.
Тематика лабораторных работ охватывает вопросы построения систем с
заданным качеством переходного процесса, синтеза оптимальных и субоптимальных алгоритмов управления и обработки информации по критерию
наивысшей точности, анализа чувствительности систем при изменении параметров входных воздействий. К каждой работе даны методические рекомендации по выполнению исследований, подробное описание программы моделирования, контрольные вопросы.
6
1. СИНТЕЗ СИСТЕМЫ СТАБИЛИЗАЦИИ СУДНА ПО КУРСУ С
ЗАДАННЫМ КАЧЕСТВОМ ПЕРЕХОДНОГО ПРОЦЕССА
Цель работы: изучение процедуры параметрического синтеза законов
управления с заданием допустимого «коридора» для переходной характеристики.
1.1.Методические указания по подготовке к работе
Как известно [1, 3, 7, 8], все реальные системы автоматики нелинейны.
Физическим звеньям свойственны явления насыщения, ограничения, люфта,
гистерезиса и т.д. Теория нелинейных систем сложна, недостаточно хорошо
развита, и, по существу, их анализ всегда является приближенным. В то же
время существует целый класс нелинейных систем, которые при определенных
допущениях можно линеаризовать, т.е. рассматривать линейные математические модели вместо нелинейных. Для таких систем нелинейные свойства не
играют существенной роли, и их исследование, как правило, не представляет
особых технических трудностей.
С другой стороны, существует обширное множество систем, для которых
нелинейные свойства являются принципиально важными и применение линейных моделей приводит к качественно неверным результатам. Причем сложность и большое разнообразие свойств таких систем не позволяют использовать какой-либо универсальный подход к их анализу и синтезу. В связи с этим
при исследовании нелинейных систем широкое распространение получили
численные методы, которые могут применяться, в том числе, и с целью их оптимизации.
Настоящая лабораторная работа посвящена синтезу регулятора для нелинейной системы управления с заданным качеством переходного процесса.
Данный подход к синтезу основывается на целенаправленном поиске таких
значений настраиваемых параметров в законе управления, чтобы траектории
7
движения динамической системы не выходили за пределы заданных областей
(«коридоров») [7]. Существо подхода заключается в минимизации некоторой
целевой функции, определяющей меру выхода исследуемой траектории за пределы «коридора». Если определить эту функцию таким образом, чтобы ее нулевое значение соответствовало касанию траектории границ допустимых областей, то разрешимость задачи синтеза будет определяться наличием неположительных значений функции.
Описанный подход к параметрическому синтезу законов управления реализован в математической среде Matlab Simulink в форме пакета прикладных
программ NCD-blockset. Это инструментальное средство позволяет задавать
границы допустимого «коридора» в визуальном режиме, обеспечивает вычисление меры выхода за его пределы, осуществляет запуск численного метода
решения оптимизационной задачи и находит ее решение, если оно существует.
Рассмотрим методику использования пакета NCD-blockset на примере
синтеза системы стабилизации морского судна по курсу (авторулевого).
Укрупненная структурная схема системы представлена на рис. 1.1, на котором
обозначены: УУ – устройство управления (регулятор), ИУ – исполнительное
устройство (рулевой привод), ОУ – объект управления (судно); Kз – заданное
значение курса, K – текущее значение курса, = Kз – K – отклонение курса от
заданного значения (ошибка регулирования), u – сигнал управления, – угол
перекладки руля.
Kз
+
УУ
u
ИУ
ОУ
K
–
Рис. 1.1. Структурная схема авторулевого
8
Для упрощенного описания движения судна (без учета внешних возмущений) будем использовать нелинейное дифференциальное уравнение второго
порядка в форме Номото [9, 10]:
(1 2 )
c2 c33 k1 k13 ,
1 2
(1.1)
где = ̇ – угловая скорость изменения курса.
В частном случае, если квадратичным и кубическим слагаемыми в уравнении (1.1) можно пренебречь, система стабилизации становится линейной с
передаточной функцией (от к )
WÎÓ ( s)
k1 (1 3 s)
.
s(1 1s)(1 2 s)
Установившееся значение угловой скорости при постоянном угле перекладки руля можно получить, приравняв нулю все производные в уравнении
(1.1):
c2 c3 3 k1.
(1.2)
Выражение (1.2) определяет зависимость (), графическое представление которой называется диаграммой поворотливости (рис. 1.2).
9
2
п
0
1
Рис. 1.2. Диаграммы поворотливости
Характер диаграммы поворотливости зависит от устойчивости судна на
курсе, определяемой знаками коэффициентов в уравнении (1.1). Когда судно
теоретически устойчиво на курсе, все коэффициенты в этом уравнении положительны (кривая 1). Если судно не обладает устойчивостью прямолинейного
движения, то 1 < 0, k1 < 0, c2 ≤ 0, а 2 > 0, 3 > 0 (кривая 2). Знак c3 зависит не
только от устойчивости судна, но и от вида диаграммы поворотливости при
больших перекладках руля [9]. Из опыта судовождения следует [11], что при
c2 ≤ 0 значения | | < |0| у судов неустойчивы и могут быть получены только
в результате некоторого процесса регулирования движения. Поэтому на
рис. 1.2 кривая на данном участке изображена штриховой линией. При этом
изменение угловой скорости на предельных углах обратной поворотливости
±п происходит не плавно, а скачком, в результате чего зависимость () приобретает вид петли гистерезиса.
Для анализа математической модели движения судна по курсу перейдем
от дифференциального уравнения (1.1) к описанию системы в пространстве состояний:
10
x1 x2 b1,
c
2
c
1
x 2 1
x2
x1 2 x1 x1 3 x13 b2,
1 2
1 2
1 2
1 2
(1.3)
x3 x1.
В данной модели угловой скорости соответствует переменная состояния x1, а углу – переменная состояния x3. Коэффициенты b1,2 определяются из
соотношений [1]:
b1
k13
,
1 2
13 2 3
b2 k1 1 2
.
12 22
Работа системы стабилизации курса считается удовлетворительной, когда
судно достаточно быстро и с малым перерегулированием выходит на заданный
курс и стабилизируется на нем с высокой точностью удержания (малым рысканием). Обеспечить выполнение требований быстродействия и малой колебательности
возможно
при
использовании
пропорционально-интегрально-
дифференциального закона управления (ПИД-регулятора) [9-11]:
k è dt ,
u kï k ä
(1.4)
где kп, kд, kи – коэффициенты регулятора. При этом управление по интегралу от
углового рассогласования обеспечивает нулевую статическую ошибку регулирования, а управление по производной служит для демпфирования собственных колебаний системы [7, 9].
Выражению (1.4) соответствует передаточная функция регулятора
11
WÓÓ ( s) kï k ä s kè s .
Поскольку использование интегральной составляющей в законе управления не способствует улучшению качества переходного процесса, ее целесообразно включать только в режиме стабилизации. Поэтому при отклонениях курса от заданного значения на величину 5 и более будем полагать kи = 0.
Управление u отрабатывается рулевым приводом, математическую модель которого примем в виде [11]
u ,
( min sign ) ( min ),
sat( , ) ( ),
4
max
(1.5)
max
где min – минимальный угол перекладки руля, отрабатываемый приводом (половина ширины зоны нечувствительности), max – максимальный угол перекладки
руля,
max
–
максимальная
скорость
перекладки
руля,
5 – постоянная времени привода; sat(x, xmax) – функция насыщения:
x, x xmax ,
sat( x, xmax ) xmax , x xmax ,
x , x x ;
max
max
(x) – единичная функция Хевисайда:
1, x 0,
( x)
0, x 0.
12
При min u max, u max система (1.5) допускает линеаризованное
представление:
(u ) 4 ,
что позволяет рассматривать привод как апериодическое звено первого порядка с передаточной функцией
WÈÓ ( s)
1
.
1 4 s
Однако в том случае, когда указанные условия не выполняются, существенную роль играют присущие рулевому приводу нелинейности типа насыщения и нечувствительности. Первое из указанных свойств привода проявляется в переходном режиме системы управления, второе – в режиме стабилизации.
Зададим «коридор» для переходной характеристики замкнутой системы
исходя из предельно допустимых значений перерегулирования m и длительности переходного процесса tп (рис. 1.3). При этом определим длительность переходного процесса как время, по истечении которого величина курса укладывается в диапазон 2% относительно заданного значения Kз.
Задача синтеза заключается в нахождении таких значений коэффициентов регулятора kп, kд, kи, чтобы переходная характеристика замкнутой системы
укладывалась в допустимый «коридор». Остальные параметры системы считаются заданными.
13
K/Kз
1+m
1,02
1
0,98
0
tп
t
Рис. 1.3. Допустимый «коридор» для переходной характеристики
1.2.Описание программы моделирования
Лабораторная работа выполняется в среде Matlab 6.5 с использованием
пакета Simulink.
В соответствии с уравнениями (1.3)-(1.5) формируется Simulink-модель,
показанная на рис. 1.4. Текущим отклонением курса от заданного значения
Kз является выход блока Integrator2, при этом входной сигнал Kз задается в виде начального условия в блоке Integrator2.
Текущий курс (вход блока Gain2) формируется как разность между курсом заданным и величиной : K = Kз – . Блок Gain2 служит для нормировки
курса к заданному значению с той целью, чтобы установившееся значение
входного сигнала блока NCD Outport было равно единице.
14
В нулевой момент времени текущее значение курса для определенности
принимается равным нулю. При том, что начальные условия для ненулевые
( = Kз), система управления будет стремиться свести к нулю значение
= Kз – K, в результате чего величина курса изменяется от нуля до заданного
значения. Тем самым осуществляется маневр по курсу на величину Kз относительно первоначального значения.
b1
K3
Gain
b2
Gain1
Constant
1
s
1
s
1
s
Integrator
Integrator1
Integrator2
MATLAB
Function
MATLAB Fcn
Ground
Gain3
-K-
Gain4
kp
Switch
Gain5
MATLAB
Function
kd
MATLAB Fcn1
1
s
Integrator4
Gain6
Saturation
Gain7
1
s
Integrator3
-K-
NCD
OutPort 1
-KDead Zone
ki
Gain2
NCD Outport
R2D
Scope
Radians
to Degrees
Рис. 1.4. Модель исследуемой системы в пакете Simulink
В блоке NCD Outport задается допустимый «коридор» для переходной
характеристики. После запуска процедуры настройки параметров регулятора в
графической области отображается ход ее выполнения в виде исходного графика переходного процесса белого цвета и текущего изменяющегося графика
зеленого цвета. В это же время в командном окне Matlab отображаются текущие значения настраиваемого параметра и целевой функции оптимизации
15
(функции штрафа). По окончании выполнения процедуры оптимальные значения настраиваемых переменных, соответствующие кривой зеленого цвета, сохраняются в рабочем пространстве Matlab.
Отметим, что поскольку в пакете NCD-Blockset используется градиентный метод поиска минимума функции штрафа, попадание в точку глобального
экстремума не гарантируется, даже если на заданных ограничениях он в принципе достижим. В этом случае следует изменить начальное приближение для
настраиваемого параметра и повторить процедуру настройки.
1.3.Порядок выполнения лабораторной работы
1. Получить вариант задания у преподавателя.
2. Запустить Matlab, создать новый mdl-файл.
3. Сформировать Simulink-модель в соответствии с рис. 1.4. Ввести:
в блок Gain
Gain:
b1
в блок Constant
Constant value:
K3
в блок Gain1
Gain:
b2
в блоки Integrator, Integrator3
Initial condition:
0
в блок Integrator1
Initial condition:
0
(если судно устойчиво на курсе);
Initial condition:
1/с2
(если судно неустойчиво на курсе и c3 = 0 либо
Initial condition:
c3 c2 );
-1/sqrt(-c3)
(если судно неустойчиво на курсе и c2 = 0)
в блок Integrator2
16
Initial condition:
K3
в блок Gain2
Gain:
1/K3
в блок MATLAB Fcn
MATLAB Function: abs(u)
в блок Gain3
Gain:
(tau1+tau2)/tau1/tau2
в блок Gain4
Gain:
kp
в блок Switch
Criteria for passing first input:
Threshold:
u2 >= Threshold
5*pi/180
в блок Gain5
Gain:
ki
в блок MATLAB Fcn1
MATLAB Function: 1/tau1/tau2*(u + c2*u*abs(u) + c3*u^3)
в блок Gain6
Gain:
kd
в блок Integrator4
Initial condition:
0
Limit output
Upper saturation limit:
dmax
Lower saturation limit:
- dmax
в блок Saturation
Upper limit:
ddmax
Lower limit:
- ddmax
в блок Gain7
Gain:
1/tau4
в блок Dead Zone
17
Start of dead zone:
-dmin
End of dead zone:
dmin
Задаваемые начальные условия соответствуют прямолинейному движению с установившимся курсом (для устойчивых судов) либо движению с постоянной скоростью циркуляции (для неустойчивых судов).
4. В меню Simulation выбрать пункт Simulation Parameters и установить время моделирования 300 с.
5. Создать m-файл в соответствии с приложением 1. Ввести значения
параметров в соответствии с вариантом задания (табл. 1.1), обращая внимание
на их размерности.
Табл. 1.1. Варианты заданий
Ваk1 ,
1, 2, 3, c2,
ри–1
с
с
с
с
с
ант
1
0,03 10 1 2
20
2
0,05 20 2 5
35
3
0,1
30 3 7
20
4
–0,1 –40 4 10
0
5 –0,12 –50 5 10 –50
6 –0,13 –60 6 15
0
7 –0,15 –70 7 20 –50
c3 , с2
min, max, ̇ max, 4,
град град град/с с
Kз, m, tп,
град % с
20
0
40
–900
10
–700
0
0,5
0,3
0,4
0,3
0,3
0,4
0,3
60
60
90
90
90
90
90
35
25
35
35
35
35
25
2
3
4
3
3
4
2
0,3
0,3
0,25
0,4
0,3
0,3
0,4
5 100
10 90
10 60
10 90
10 90
20 75
20 90
6. Вернуться в окно созданной Simulink-модели. Двойным щелчком
мыши по блоку NCD Outport раскрыть графическое окно, представленное на
рис. 1.5.
18
Рис. 1.5. Графическое окно блока NCD Outport
7. Войти в меню Options, выбрать пункт Time Range… и установить
значение tmax = 300с.
8. Построить «коридор», в пределах которого должен находится выходной сигнал исследуемой системы, в соответствии с заданными величинами перерегулирования m и длительности переходного процесса tп. Это можно сделать, передвигая красные линии, являющиеся границами «коридора», при помощи мыши. Кроме того, местоположение этих линий можно установить точно
(не в визуальном режиме) при помощи диалоговой панели Constraint Editor,
возникающей при щелчке правой кнопкой мыши по красной линии. В качестве
примера на рис. 1.6 приведен вид допустимого «коридора», соответствующий
перерегулированию 10% и длительности переходного процесса 60 с. Положение границ при этом задается в виде [0 -0.01 60 -0.01]; [0 1.1 60 1.1];
[60 0.98 300 0.98]; [60 1.02 300 1.02]. Небольшое отрицательное значение координаты местоположения самой нижней из границ «коридора» задается для
«облегчения» процедуры настройки, так как в противном случае (при установлении нижней границы точно на нулевом уровне) в момент времени, равный
нулю, переходная характеристика неизбежно касается нижней границы «кори-
19
дора», в связи с чем функция штрафа лишается возможности принимать отрицательные значения.
Рис. 1.6. Построение допустимого «коридора» для переходной
характеристики
9. В окне созданного mdl-файла войти в меню Simulation Parameters и
установить время моделирования, равное tmax.
10. В блоке NCD Outport войти в меню Optimization, выбрать пункт Parameters…, ввести:
Tunable Variables:
kp kd ki
Lower bounds:
-100 -100 -0.1
Upper bounds:
0 0 -0.001
Discretization Interval: 1
Variable Tolerance:
0.01
Constraint Tolerance:
0.01
Display optimization information
Stop optimization as soon as the constraints are achieved
Compute gradients with better accuracy (slower)
Зафиксировать введенные значение нажатием кнопки Done.
20
11. Запустить процедуру настройки коэффициентов регулятора, нажав
виртуальную кнопку Start.
12. Наблюдать развитие процесса настройки коэффициентов по изменению кривой переходной характеристики (рис. 1.7). Одновременно контролировать ход изменения значения функции штрафа MAX{g} (максимального выхода кривой переходного процесса за границы «коридора»), отображаемой в командном окне Matlab.
Рис. 1.7. Отображение хода процесса определения настраиваемых переменных
13. Убедиться в завершении процедуры настройки по появлению сообщения "Optimization Converged Successfully" в командном окне Matlab.
14. Ввести в командном окне Matlab
>> [kp ki kd]
и зафиксировать полученные значения настраиваемых переменных.
15. В случае если значение MAX{g}, соответствующее последней итерации, оказалось отрицательным (или равным нулю), считать задачу синтеза разрешенной, а полученные значения настраиваемых переменных – значениями
коэффициентов регулятора kп, kи, kд. В противном случае откорректировать
начальные приближения (как правило, следует на 10-20% увеличить по модулю
21
величину kд по отношению к полученному значению), введя новое приближение в командной строке Matlab (kd = …), и повторить процедуру синтеза,
начиная с п. 11.
16. В окне mdl-файла нажать кнопку Start Simulation. Двойным щелчком
мыши по блоку Scope раскрыть окно со строящимся графиком переходного
процесса. Изменив масштаб графика, определить амплитуду рыскания в режиме стабилизации.
1.4.Оформление отчета
Отчет о лабораторных исследованиях должен содержать:
структурную схему исследуемой системы;
модель системы в пакете Simulink;
исходные данные для синтеза;
полученные значения коэффициентов регулятора и соответствующий
им график переходного процесса;
полученное значение амплитуды рыскания в режиме стабилизации.
1.5.Контрольные вопросы
1. Поясните суть процедуры определения настраиваемых переменных с
использованием пакета NCD-Blockset.
2. Поясните процедуру линеаризации модели объекта управления.
3. Поясните процедуру линеаризации модели рулевого привода.
4. Нелинейности какого вида характерны для рулевого привода? Каким
образом они проявляются в рассматриваемой системе?
5. Запишите передаточную функцию замкнутой линеаризованной системы.
6. Каким образом обеспечивается астатизм в рассматриваемой системе?
7. Является ли используемый в системе регулятор линейным? Почему?
22
8. Чему равно установившееся значение входного сигнала блока NCD
Outport? Почему?
9. Как по графику переходной характеристики определить перерегулирование и длительность переходного процесса?
10.В схеме на рис. 1.4 покажите точки, где можно наблюдать сигналы
x1(t), x2(t), x3(t), u(t), (t).
23
2. ИССЛЕДОВАНИЕ ТОЧНОСТИ ОПРЕДЕЛЕНИЯ УГЛОВ
ОРИЕНТАЦИИ ПОДВИЖНОГО ОБЪЕКТА С ИСПОЛЬЗОВАНИЕМ
БЕСПЛАТФОРМЕННОЙ ИНЕРЦИАЛЬНОЙ ВЕРТИКАЛИ
Цель работы: ознакомление с алгоритмом работы бесплатформенной
инерциальной вертикали.
2.1.Методические указания по подготовке к работе
Бесплатформенная инерциальная вертикаль предназначена для выработки углов ориентации маневренных объектов (МО). Основой вертикали является
измерительный блок, содержащий три линейных акселерометра, измеряющих
составляющие кажущегося ускорения, и три датчика угловой скорости (ДУС),
измеряющих составляющие угловой скорости вращения. Выходные сигналы
датчиков поступают непосредственно в устройство оценивания (УО), которое
определяет мгновенное направление осей чувствительности акселерометров в
опорной системе координат. Укрупненная структурная схема системы представлена на рис. 2.1. Штриховыми линиями на схеме обозначены связи, которые далее будут искусственно введены для получения математической модели,
удобной при синтезе устройства оценивания.
v
w, , J, g
x4, x5, x6
МО
Акс
ДУС
wп
п
УО
˄ g˄
˄ J,
,
x1, x2, x3
Рис. 2.1. Структурная схема инерциальной вертикали
24
При построении вычислительного алгоритма бесплатформенной системы
ориентации примем за основу уравнения Пуассона [12, 13], описывающие изменение матрицы направляющих косинусов
cos cos J
sin J
sin cos J
C sin sin g cos sin J cos g cos J cos g cos sin g sin sin J cos g ,
sin cos g cos sin J sin g cos J sin g cos cos g sin sin J sin g
где J – угол тангажа, – угол рыскания, g – угол крена.
В матричной форме эти уравнения имеют вид
dCT
CT ω,
dt
(2.1)
где – кососимметрическая матрица
0
ω z
y
z
0
x
y
x ,
0
составленная из проекций вектора угловой скорости объекта на оси связанной с
ним системы координат.
Выражение (2.1) в скалярной форме представляет собой следующую систему дифференциальных уравнений:
C11 C21 z C31 y ,
(2.2)
C 21 C31 x C11 z ,
(2.3)
C 31 C11 y C21 x ,
(2.4)
25
C12 C22 z C32 y ,
(2.5)
C 22 C32 x C12 z ,
(2.6)
C 32 C12 y C22 x ,
(2.7)
C13 C23 z C33 y ,
(2.8)
C 23 C33 x C13 z ,
(2.9)
C 33 C13 y C23 x .
(2.10)
Очевидно, что каждая тройка уравнений Пуассона (2.2)-(2.4), (2.5)-(2.7) и
(2.8)-(2.10) интегрируется независимо от остальных. При этом нахождение девяти направляющих косинусов даст избыточную информацию для определения
параметров ориентации, поэтому на практике в бортовых вычислителях интегрируют не все девять уравнений, а только ту их часть, которая позволяет
наиболее простым путем вычислить углы ориентации. Далее ограничимся рассмотрением уравнений (2.5)-(2.10) с начальными условиями
C12 (0) sin J0 ,
C22 (0) cos J0 cos g 0 ,
C32 (0) cos J0 sin g 0 ,
C13 (0) sin 0 cos J0 ,
C23 (0) cos 0 sin g 0 sin 0 sin J0 cos g 0 ,
C33 (0) cos 0 cos g 0 sin 0 sin J0 sin g 0 .
Здесь J0, 0, g0 – начальные углы ориентации.
Недостающие элементы матрицы направляющих косинусов можно определить из соотношений
C11 C22 C33 C23C32 ,
(2.11)
26
C21 C13C32 C12 C33 ,
(2.12)
C31 C12 C23 C13C22 .
(2.13)
Тогда параметры ориентации вычисляются по формулам
J arctg
C12
2
2
C 22
C32
, arctg
C13
C
, g arctg 32 .
C11
C22
Из формул (2.2)-(2.10) следует, что для вычисления параметров ориентации необходимо знание составляющих угловой скорости вращения объекта.
Такая информация поступает с выходов ДУС с точностью до ошибок измерений:
ï
äð
ï
äð
x ïx äð
x , y y y , z z z ,
где пx , пy , пz – измеренные (приборные) значения составляющих угловой скорости вращения, являющиеся входами исследуемой динамической системы;
др
др
др
x , y , z – угловые скорости дрейфа ДУС. Методические ошибки ДУС бу-
дем считать здесь пренебрежимо малыми.
Таким образом, система уравнений относительно направляющих косинусов Cij примет следующий вид:
ï
äð
C12 C 22 ïz C22 äð
z C32 y C32 y ,
(2.14)
ï
äð
C 22 C32 ïx C32 äð
x C12 z C12 z ,
(2.15)
ï
äð
C 32 C12 ïy C12 äð
y C 22 x C 22 x ,
(2.16)
ï
äð
C13 C 23 ïz C 23 äð
z C33 y C33 y ,
(2.17)
ï
äð
C 23 C33 ïx C33 äð
x C13 z C13 z ,
(2.18)
27
ï
äð
C 33 C13 ïy C13 äð
y C 23 x C 23 x .
(2.19)
др
др
Угловые скорости дрейфа др
x , y , z будем считать центрированными
экспоненциально коррелированными случайными процессами со среднеквадратическими отклонениями 1, 2, 3 и постоянными времени корреляции T1,
T2, T3 соответственно. Уравнения формирующих фильтров для таких процессов
имеют вид [3, 14]
1 äð
1
äð
x T1 x 2T1 1x1 ,
(2.20)
1 äð
1
äð
y T2 y 2T2 2 x 2 ,
(2.21)
1 äð
1
äð
z T3 z 2T3 3 x 3 ,
(2.22)
где xi – некоррелированные друг с другом центрированные белые шумы единичной интенсивности.
Составляющие абсолютного ускорения объекта wgx, wgy, wgz также представим в виде центрированных экспоненциально коррелированных случайных
процессов со среднеквадратическими отклонениями gx, gy, gz и постоянными
времени корреляции Tgx, Tgy, Tgz соответственно:
w gx Tgx1wgx 2Tgx1 gx x 4 ,
(2.23)
w gy Tgy1wgy 2Tgy1 gy x5 ,
(2.24)
w gz Tgz1wgz 2Tgz1 gz x 6 .
(2.25)
Выходные сигналы акселерометров wпx , wпy , wпz будут представлять собой
зашумленные измерения кажущихся ускорений в проекциях на их оси чувствительности. При движении объекта с постоянной скоростью
28
где
wxï wgx C11 ( wgy g )C12 wgz C13 1 ,
(2.26)
w ïy wgx C 21 ( wgy g )C 22 wgz C 23 2 ,
(2.27)
wzï wgx C31 ( wgy g )C32 wgz C33 3 ,
(2.28)
j – некоррелированные друг с другом и с xi центрированные белые шумы
с заданными интенсивностями Rj; g – ускорение силы тяжести.
Дифференциальные уравнения (2.14)-(2.25) и алгебраические уравнения
(2.26)-(2.28) представляют собой модель динамической системы, которая может быть записана в векторно-матричной форме:
x f (x, u) Bξ,
(2.29)
z h(x) v,
(2.30)
др
др
T
где x = (C12 C22 C32 C13 C23 C33 др
x y z wgx wgy wgz) – вектор состояния;
z = (wпx
wпy
wпz )T – вектор измерений; u = (пx
пy
пz )T – вектор входных
сигналов; x = (x1 x2 x3 x4 x5 x6)T – вектор возмущений; v = (
1
2
3)
T
–
вектор шумов измерений; B – матрица возмущений размерности 123, ненулевые элементы которой равны
B7,1 2T11 1 , B8, 2 2T21 2 , B9,3 2T31 3 ,
B10, 4 2Tgx1 gx , B11,5 2Tgy1 gy , B12,6 2Tgz1 gz .
Компоненты вектор-функции f(x,u) будут иметь вид
f1 (x, u) x2 x9 x3 x8 x2u3 x3u2 ,
f 2 (x, u) x3 x7 x1 x9 x3u1 x1u3 ,
29
f 3 (x, u) x1 x8 x2 x7 x1u2 x2u1 ,
f 4 (x, u) x5 x9 x6 x8 x5u3 x6u2 ,
f 5 (x, u) x6 x7 x4 x9 x6u1 x4u3 ,
f 6 (x, u) x4 x8 x5 x7 x4u2 x5u1 ,
f 7 (x, u) T11 x7 ,
f 8 (x, u) T21 x8 ,
f 9 (x, u) T31 x9 ,
f10 (x, u) Tgx1 x10 ,
f11 (x, u) Tgy1 x11 ,
f12 (x, u) Tgz1 x12 .
Компоненты вектор-функции h(x) будут иметь вид
h1 (x) x10 ( x2 x6 x3 x5 ) ( x11 g ) x1 x12 x4 ,
h2 (x) x10 ( x3 x4 x1 x6 ) ( x11 g ) x2 x12 x5 ,
h3 (x) x10 ( x1 x5 x2 x4 ) ( x11 g ) x3 x12 x6 .
Очевидно, функции f(x,u), h(x) являются нелинейными, т.к. их компоненты содержат произведения переменных состояния. Поэтому для решения задачи фильтрации (оценивания компонент вектора состояния) будем использовать
обобщенный нелинейный фильтр Калмана-Бьюси. Уравнения фильтра записываются в виде [14, 15]
xˆ f (xˆ , u) K (z h(xˆ )),
(2.31)
f (xˆ , u) 1
R ,
x T
(2.32)
KP
30
T
T
P f (xˆ , u) P P f (xˆ , u) P h (xˆ ) R 1 h(xˆ ) P BB T ,
x
x
x T
x T
(2.33)
где x̂ – текущая оценка вектора состояния; P – матрица ковариаций вектора
ошибок ε x x̂; R – матрица интенсивностей шумов измерений (диагональная
матрица, составленная из элементов Rj); K – матричный коэффициент усиления
фильтра.
Здесь
f (xˆ , u)
x
T
представляет собой матрицу вида
f1
x1
f
f (x, u) 2
x1
x T
f12
x
1
вычисленную в точке x x̂, а
f1
x2
f 2
x2
...
f12
x2
f1
x12
f 2
...
x12 ,
f12
...
x12
...
f T (xˆ , u)
– транспонированную матрицу. Аналоx
h(xˆ ) h T (xˆ )
.
гично определяются производные
,
x
x T
В данном случае
f (x, u)
x T
31
0
x9 u3 x8 u 2 0
0
0
0 x3 x2 0 0 0
0
x7 u1 0
0
0
x3 0 x1 0 0 0
x9 u3
x u x u
0
0
0
0 x2 x1 0 0 0 0
2
7
1
8
0
0
0
0
x9 u3 x8 u 2 0 x6 x5 0 0 0
0
0
0 x9 u 3
0
x7 u1 x6 0 x4 0 0 0
0
0
0 x8 u 2 x7 u1
0 x5 x4 0 0 0 0
0
0
0
0
0
0 T11 0 0 0 0 0 ,
0
0
0
0
0
0
0 T21 0 0 0 0
0
0
0
0
0
0
0
0 T31 0 0 0
1
0
0
0
0
0
0
0
0 0 Tgx 0 0
0
0
0
0
0
0
0
0 0 0 Tgy1 0
1
0
0
0
0
0
0
0
0 0 0 0 Tgz
h(x)
x T
x11 g
x6 x10
x x
5 10
x6 x10
x11 g
x4 x10
x5 x10
x4 x10
x11 g
x12
x3 x10
x2 x10
x3 x10
x12
x1 x10
x2 x10
x1 x10
x12
0 0 0 x 2 x6
0 0 0 x3 x 4
0 0 0 x1 x5
x1
x2
x3
x4
x5 .
x6
Начальными условиями для интегрирования уравнений (2.31)-(2.33) являются априорная оценка вектора состояния x̂(0) и априорная матрица ковариаций вектора ошибок P(0) = P0, где P0 – заданная априорная матрица ковариаций вектора состояния (диагональная матрица, состоящая из дисперсий компо2
нент вектора состояния 0i
). Оптимальной априорной оценкой вектора состоя-
ния является его безусловное математическое ожидание, которое при малых
углах
J0,
0,
g0
будем
считать
приближенно
равным
вектору
xˆ (0) (0 1 0 0 0 1 0 0 0 0 0 0 0) T . Диагональные элементы априорной мат-
рицы ковариаций P0 определим как
32
2
2
2
2
2
01
02
03
04
05
10 4 ;
2
2
2
2
06
10 6 ; 07
12 ; 08
22 ; 09
32 ;
02 10 2gx ; 02 11 2gy ; 02 12 2gz .
2.2. Описание программы моделирования
Лабораторная работа выполняется в среде Matlab 6.5 с использованием
пакета Simulink.
Для анализа точности оценивания углов ориентации формируется Simulink-модель в соответствии с рис. 2.2. Модели блоков Object (объект) и Kalman
Filter (фильтр Калмана), сформированные как подсистемы (блоки типа Subsystem) в соответствии с уравнениями (2.14)-(2.25) и (2.31)-(2.33), показаны на
рис. 2.3 и 2.4 соответственно.
Исходные данные и начальные условия для интегрирования систем дифференциальных уравнений размещаются в m-файле.
Возмущающие входные воздействия моделируются широкополосными
процессами с малой постоянной времени (0,01 с).
Движение объекта задается посредством формирования составляющих
угловой скорости, изменяющихся по пилообразному закону с заданными амплитудами x, y, z, частотами x, y, z и случайной начальной фазой.
В
работе
для
всех
вариантов
приняты
значения
констант
R1 = R2 = R3 = 10–7 м2/с3; T1 = T2 = T3 = 10 с; Tgx = Tgy = Tgz = 1 с; период изменения
составляющих
угловой
скорости
движения
объекта
10 с
(x = y = z = /5 с–1).
33
x
MATLAB
Demux
Function
u
angles
e(C12),+/-3sigma
z
Object
Demux
e(C22),+/-3sigma
e(C32),+/-3sigma
u
x^
z
P
e(C13),+/-3sigma
Kalman Filter
e(C23),+/-3sigma
MATLAB
Function
MATLAB
Function
e(C33),+/-3sigma
sigma
angles est
Scope1
e(omegdx),+/-3sigma
e(omegdy ),+/-3sigma
Demux
theta, deg
e(omegdz),+/-3sigma
e(wgx),+/-3sigma
Demux
psi, deg
e(wgy ),+/-3sigma
gamma, deg
e(wgz),+/-3sigma
Scope
Scope2
Рис. 2.2. Главное окно модели
34
K*u
Band-Limited
White Noise
Gain3
1
s
Sine Wave
Sign
-K-
Integrator
1
s
Sine Wave1
Sign1
1
s
Sine Wave2
Sign2
Gain
-K-
Integrator1
omega_y
Gain1
-K-
Integrator2
omega_x
omega_z
Gain2
2
MATLAB
Function
u(t)
u(t)
u
MATLAB
Function
MATLAB Fcn
f(x,u)
x(t)
dx(t)/dt
1
s
Integrator3
MATLAB
Function
h(x)
1
3
z
x
Band-Limited
White Noise1
Рис. 2.3. Блок Object
35
Reshape
1
Matrix
Multiply
Product
1
s
u
MATLAB
Function
x^
df(x^,u)/dx'
Matrix
Multiply
Reshape
Integrator
2
MATLAB
Function
P
1
s
f(x^,u)
P
P
Integrator1
u
1
x^
2
z
MATLAB
Function
MATLAB
Function
h(x^)
dh(x^)/dx'
Product2
T
Matrix
Multiply
Math
Product3
Function
T
Reshape
u
Reshape1
Math
Function1
Matrix
Multiply
R^-1
Constant
Product1
B*B'
Constant1
dP/dt
Рис. 2.4. Блок Kalman Filter
2.3.Порядок выполнения лабораторной работы
1. Получить вариант задания у преподавателя.
2. Запустить Matlab, создать новый mdl-файл.
3. Сформировать Simulink-модель в соответствии с рис. 2.2. Ввести:
в блок angles
Matlab Function:
[atan(u(1)/sqrt(u(2)^2+u(3)^2));
-atan(u(4)/(u(2)*u(6)-u(3)*u(5))); -atan(u(3)/u(2))]*180/pi
Output dimensions: 3
в блок angles est
Matlab Function:
[atan(u(1)/sqrt(u(2)^2+u(3)^2));
-atan(u(4)/(u(2)*u(6)-u(3)*u(5))); -atan(u(3)/u(2))]*180/pi
Output dimensions: 3
в блок sigma
Matlab Function:
[-3*sqrt(diag(u));3*sqrt(diag(u))]
36
Output dimensions: 24
в блок Bus Creator
Number of inputs:
36
в блоки Bus Selector (сверху вниз):
Selected signals:
signal1
signal13
signal25
Selected signals:
signal2
signal14
signal26
Selected signals:
signal3
signal15
signal27
Selected signals:
signal4
signal16
signal28
Selected signals:
signal5
signal17
signal29
Selected signals:
signal6
signal18
signal30
Selected signals:
signal7
signal19
signal31
Selected signals:
signal8
signal20
signal32
Selected signals:
signal9
signal21
signal33
Selected signals:
signal10
signal22
signal34
Selected signals:
signal11
signal23
signal35
Selected signals:
signal12
signal24
signal36
Muxed output (во всех блоках Bus Selector)
в блок Scope (меню ′Scope′ parameters)
вкладка General:
Number of axes: 3
Time range:
50
вкладка Data history:
Limit data points to last:
в блоки Scope1, Scope2 (меню ′Scope′ parameters)
вкладка General:
Number of axes: 6
Time range:
50
вкладка Data history:
Limit data points to last:
37
В блоках Object и Kalman Filter задать необходимое количество входов и
выходов в соответствии с рис. 2.2.
4. Раскрыть блок Object. Сформировать Simulink-модель в соответствии с
рис. 2.3. Ввести:
в блок Band-Limited White Noise
Noise power: 1
Sample time: 0.01
в блок Sine Wave
Amplitude:
Bias:
1
0
Frequency (rad/sec): Omegx
Phase (rad):
fi
Sample time: 0
в блок Sine Wave1
Amplitude:
Bias:
1
0
Frequency (rad/sec): Omegy
Phase (rad):
fi
Sample time: 0
в блок Sine Wave2
Amplitude:
Bias:
1
0
Frequency (rad/sec): Omegz
Phase (rad):
fi
Sample time: 0
в блок Integrator
Initial condition:
(-0.5+abs(fi/pi))*pi/Omegx
в блок Integrator1
Initial condition:
(-0.5+abs(fi/pi))*pi/Omegy
38
в блок Integrator2
Initial condition:
(-0.5+abs(fi/pi))*pi/Omegz
в блок Gain
Gain:
Omeg_x*2*Omegx/pi
в блок Gain1
Gain:
Omeg_y*2*Omegy/pi
в блок Gain2
Gain:
Omeg_z*2*Omegz/pi
в блок Gain3
Gain:
B
Multiplication:
Matrix(K*u)
в блок MATLAB Fcn
Matlab Function:
[u(1)-u(10); u(2)-u(11); u(3)-u(12)]
Output dimensions: 3
в блок f(x,u)
Matlab Function:
[u(5)*u(12)-u(6)*u(11)+u(5)*u(3)-u(6)*u(2);
u(6)*u(10)-u(4)*u(12)+u(6)*u(1)-u(4)*u(3);
u(4)*u(11)-u(5)*u(10)+u(4)*u(2)-u(5)*u(1);
u(8)*u(12)-u(9)*u(11)+u(8)*u(3)-u(9)*u(2);
u(9)*u(10)-u(7)*u(12)+u(9)*u(1)-u(7)*u(3);
u(7)*u(11)-u(8)*u(10)+u(7)*u(2)-u(8)*u(1); -1/T1*u(10);
-1/T2*u(11); -1/T3*u(12); -1/Tgx*u(13); -1/Tgy*u(14); -1/Tgz*u(15)]
Output dimensions: 12
в блок h(x)
Matlab Function:
[u(10)*(u(2)*u(6)-u(3)*u(5))+(u(11)+g)*u(1)+u(12)*u(4);
u(10)*(u(3)*u(4)-u(1)*u(6))+(u(11)+g)*u(2)+u(12)*u(5);
u(10)*(u(1)*u(5)-u(2)*u(4))+(u(11)+g)*u(3)+u(12)*u(6)]
Output dimensions: 3
39
в блок Integrator3
Initial condition:
x0
в блок Band-Limited White Noise
Noise power: diag(R)
Sample time: 0.01
Seed:
[1111 2222 3333]
(вектор из трех произвольных целых чисел)
5. Раскрыть блок Kalman Filter. Сформировать Simulink-модель в соответствии с рис. 2.4. Ввести:
в блок Constant
Constant value:
R^-1
в блоки Product, Product1, Product2, Product3
Multiplication:
Matrix(*)
в блок Integrator
Initial condition:
xo0
в блок f(x^,u)
MATLAB Function: [u(5)*u(12)-u(6)*u(11)+u(5)*u(3)-u(6)*u(2);
u(6)*u(10)-u(4)*u(12)+u(6)*u(1)-u(4)*u(3);
u(4)*u(11)-u(5)*u(10)+u(4)*u(2)-u(5)*u(1);
u(8)*u(12)-u(9)*u(11)+u(8)*u(3)-u(9)*u(2);
u(9)*u(10)-u(7)*u(12)+u(9)*u(1)-u(7)*u(3);
u(7)*u(11)-u(8)*u(10)+u(7)*u(2)-u(8)*u(1); -1/T1*u(10);
-1/T2*u(11); -1/T3*u(12); -1/Tgx*u(13); -1/Tgy*u(14); -1/Tgz*u(15)]
Output dimensions: 12
в блок h(x^)
MATLAB Function:
[u(10)*(u(2)*u(6)-u(3)*u(5))+(u(11)+g)*u(1)+u(12)*u(4);
u(10)*(u(3)*u(4)-u(1)*u(6))+(u(11)+g)*u(2)+u(12)*u(5);
u(10)*(u(1)*u(5)-u(2)*u(4))+(u(11)+g)*u(3)+u(12)*u(6)]
40
Output dimensions: 3
в блок df(x^,u)/dx'
MATLAB Function: [0; -u(12)-u(3); u(11)+u(2); 0; 0; 0; 0; 0; 0; 0; 0; 0;
u(12)+u(3); 0; -u(10)-u(1); 0; 0; 0; 0; 0; 0; 0; 0; 0; -u(11)-u(2); u(10)+u(1);
0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; -u(12)-u(3); u(11)+u(2); 0; 0; 0; 0; 0;
0; 0; 0; 0; u(12)+u(3); 0; -u(10)-u(1); 0; 0; 0; 0; 0; 0; 0; 0; 0; -u(11)-u(2);
u(10)+u(1); 0; 0; 0; 0; 0; 0; 0; 0; u(6) ; -u(5); 0; u(9) ; -u(8); -1/T1; 0; 0; 0;
0; 0; -u(6); 0; u(4); -u(9); 0; u(7); 0; -1/T2; 0; 0; 0; 0; u(5); -u(4); 0; u(8);
-u(7); 0; 0; 0; -1/T3; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; -1/Tgx; 0; 0; 0; 0; 0; 0;
0; 0; 0; 0; 0; 0; -1/Tgy; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; -1/Tgz]
Output dimensions: 144
в блок Integrator1
Initial condition:
P0
в блок dh(x^)/dx'
MATLAB Function: [u(11)+g; -u(6)*u(10); u(5)*u(10); u(6)*u(10);
u(11)+g; -u(4)*u(10); -u(5)*u(10); u(4)*u(10); u(11)+g; u(12); u(3)*u(10);
-u(2)*u(10); -u(3)*u(10); u(12); u(1)*u(10); u(2)*u(10); -u(1)*u(10);
u(12); 0; 0; 0; 0; 0; 0; 0; 0; 0; u(2)*u(6); u(3)*u(4); u(1)*u(5); u(1); u(2);
u(3); u(4); u(5); u(6)]
Output dimensions: 36
в блок Reshape
Output dimensionality:
Customize
Output dimensions: [12,12]
в блок Reshape1
Output dimensionality:
Customize
Output dimensions: [3,12]
в блоки Math Function, Math Function1
Function:
transpose
в блок Constant1
41
Constant value:
B*B'
6. В меню Simulation выбрать пункт Simulation Parameters и установить
время моделирования 50 с.
7. Создать m-файл в соответствии с приложением 2. Ввести значения параметров в соответствии с вариантом задания (табл. 2.1), обращая внимание на
их размерности.
Табл. 2.1. Варианты заданий
Ва1, 2, 3, gx, gy, gz,
риград/с
м/с2
ант
1
0,01
0,001
2
0,02
0,002
3
0,01
0,005
4
0,005
0,001
5
0,005
0,002
6
0,005
0,005
7
0,02
0,001
x,
град/с
y,
град/с
z,
град/с
6
5
5
8
8
10
8
8
7
10
6
5
15
10
10
12
14
16
18
20
10
J0,
0,
g0,
град град град
1
0
–1
–2
1
–1
–1
–1
1
–1
1
–1
1
0
1
–1
–1
0
–1
–1
1
8. Нажатием клавиши [F5] запустить m-файл на выполнение.
9. Вернуться в главное окно Simulink-модели. Нажатием кнопки Start
Simulation запустить mdl-файл на выполнение.
10. Двойным щелчком мыши по блокам Scope1, Scope2, Scope3 раскрыть
окна со строящимися графиками:
углов ориентации J, , g и их оценок;
ошибок оценок переменных состояния с оценками точности (3).
11. Сравнить оценки углов ориентации с их истинными значениями. Записать величины максимальных ошибок для отчета.
12. Убедиться в соответствии ошибок оценок переменных состояния
оценкам их точности (проконтролировать попадание значений ошибок в коридор 3). Сделать вывод о том, каким образом ошибки зависят от ориентации
объекта.
42
13. В блоках Band-Limited White Noise, Band-Limited White Noise1 установить интенсивности входных шумов, равные нулю. Повторить процедуру
моделирования. Определить переменные состояния, переходные процессы для
которых имеют затухающий характер. По длительности переходных процессов
оценить время готовности бесплатформенной вертикали к работе.
2.4.Оформление отчета
Отчет о лабораторных исследованиях должен содержать:
структурную схему исследуемой системы;
систему дифференциальных уравнений, описывающих поведение
объекта;
полученные значения экспериментальных данных;
выводы по работе.
2.5.Контрольные вопросы
1. Какие датчики включает в себя измерительный блок бесплатформенной вертикали? Какие физические величины они измеряют?
2. Какой вид имеют корреляционные функции составляющих угловой
скорости дрейфа ДУС?
3. Какие составляющие содержат выходные сигналы акселерометров?
4. Сколько элементов матрицы направляющих косинусов дают полную
информацию об углах ориентации объекта? Какой порядок имеет исследуемая
динамическая система?
5. Какие сигналы являются входными и какие – выходными по отношению к устройству оценивания?
6. Какое распределение имеет начальная фаза в законе изменения угловой скорости движения объекта?
7. В каком из блоков модели Simulink определяются начальные условия
для интегрирования дисперсионного уравнения?
43
8. В каких блоках модели Simulink осуществляется формирование закона
изменения составляющих угловой скорости?
9. Каким образом моделируются входные возмущающие воздействия типа белого шума?
10. Является ли исследуемая система полностью наблюдаемой? Почему?
44
3. СИНТЕЗ РЕГУЛЯТОРА ДЛЯ СИСТЕМЫ УПРАВЛЕНИЯ УПРУГИМ
ОБЪЕКТОМ
Цель работы: изучение вопросов построения линейных регуляторов с
демпфирующими свойствами.
3.1.Методические указания по подготовке к работе
Стремление создать максимально легкие управляемые конструкции приводит к проявлению их упругих свойств. Учет этих свойств принципиально
важен при проектировании систем автоматического управления подвижными
объектами, подверженными значительным динамическим нагрузкам и силам
сопротивления окружающей среды. Наличие упругости способствует возникновению колебаний в системе управления на различных резонансных частотах,
в ряде случаев приводящих к потере устойчивости и разрушению конструкции.
Созданию эффективных регуляторов, подавляющих колебательные составляющие движения, препятствуют сложность получения достоверной информации об упругих свойствах нелинейного объекта и значительное изменение резонансных частот во времени. Применяющийся при решении таких задач
метод конечных элементов оказывается непригодным для аэрокосмических аппаратов сложной формы, состоящих из сотен и тысяч деталей, даже с учетом
широких возможностей современной вычислительной техники. Одним из возможных путей преодоления возникающих проблем является переход от исходной системы уравнений в частных производных, описывающих движение
упругого объекта, к линеаризованной системе уравнений большей размерности
с дальнейшим ее упрощением [16]. Далее будем считать, что математическая
модель исследуемого объекта для отдельного режима движения получена с использованием данной методики и записана в виде передаточной функции
45
W ( s)
k ( s 1 )...(s )(s 2 211s 12 )...(s 2 2 p p s 2p )
( s 1 )...(s )(s
2
2d11s 12 )...(s 2 2d s 2 )
(3.1)
с постоянными параметрами k, i, i, i, i, i, di. При этом коэффициенты i задают апериодические составляющие движения объекта, а коэффициенты i, di –
колебательные составляющие, или моды упругих колебаний. Для реальных
объектов порядок числителя передаточной функции меньше порядка знаменателя, т.е. + 2p < + 2.
Задача управления заключается в подавлении (демпфировании) мод
упругих колебаний и получении замкнутой системы с приемлемыми динамическими характеристиками. Структурная схема синтезируемой системы управления представлена на рис. 3.1, где g(t) – задающее воздействие; u(t) – управляющее воздействие; w(t) – внешнее возмущение, приведенное к входу системы;
z (t) – выходной сигнал объекта управления; z(t) – наблюдаемый (измеряемый) сигнал;
(t) – шум измерений; Wр(s) – передаточная функция регулятора.
Возмущения w(t) и
(t) будем считать некоррелированными центрированными
белыми шумами с известными интенсивностями N1 и N2.
w(t)
~
g(t)
z(t)
W(s)
u(t)
Wp(s)
z(t)
u(t)
Рис. 3.1. Структурная схема замкнутой системы
Синтез закона управления сложными динамическими системами удобно
производить с использованием аппарата пространства состояний. С целью перехода от передаточной функции объекта к его описанию в пространстве со46
стояний, в котором каждая из переменных состояния отвечала бы за ту или
иную колебательную либо апериодическую составляющую движения, предварительно представим передаточную функцию (3.1) в виде суммы простейших
дробей. При отсутствии кратных корней у знаменателя такое представление
имеет вид
W ( s)
r
r1
r
2 ... n .
s p1 s p2
s pn
Число слагаемых n = + 2 определяет порядок системы, или размерность вектора состояния.
В результате разложения передаточной функции W(s) на элементарные
дроби образуются две группы слагаемых: с вещественными и с комплексно сопряженными полюсами [6]. Рассмотрим каждую из них.
Пусть Wi ( s)
ri
– элементарная передаточная функция с вещественs pi
ным полюсом pi. Такой передаточной функции (от входа u к выходу z
(i)
) соот-
ветствует описание в пространстве состояний
xi pi xi u ,
~
z (i ) r x .
i i
Здесь z
(i)
– i-я аддитивная составляющая выходного сигнала z .
Теперь предположим, что функции Wi ( s)
ri
r
и Wi ( s) i 1 имеs pi
s pi 1
ют комплексно сопряженные полюса pi и pi+1 (соответственно, комплексно сопряженными будут также коэффициенты ri и ri+1). Чтобы избавиться от мнимой
единицы в рассматриваемых передаточных функциях, объединим эти слагаемые в одно:
47
~
Wi ( s) Wi ( s) Wi 1 ( s)
ri
r
r p r p s(ri ri 1 )
i 1 2i i 1 i 1 i
.
s pi s pi 1
s ( pi pi 1 ) s pi pi 1
Передаточной функции W
i (s)
соответствует описание в пространстве
состояний [1]
xi xi 1 ,
xi 1 pi pi 1 xi ( pi pi 1 ) xi 1 u,
~
z (i ) (ri pi 1 ri 1 pi ) xi (ri ri 1 ) xi 1.
Преобразованная модель объекта управления будет представлять собой
параллельное соединение звеньев с передаточными функциями W
(j)
(s) в соот-
ветствии с рис. 3.2. При этом каждое из таких звеньев будет формировать либо
апериодическую, либо колебательную составляющую z
(i)
выходного сигнала
z , причем в первом случае передаточной функции W
(j)
(s) будет соответ-
ствовать одна из функций Wi(s) с вещественным полюсом, а во втором – одна
i (s).
~
W(1)(s)
~
g+u+w
W(2)(s)
z~
...
из функций W
~
W(+)(s)
Рис. 3.2. Представление структуры объекта в виде
параллельного соединения звеньев
48
Таким образом, математическую модель свободного движения упругого
объекта (при отсутствии задающего воздействия) с учетом действующего на
входе системы возмущения w(t) и шума измерений
(t) можно представить в
векторно-матричной форме следующим образом:
x Ax B(u w),
z Cx υ .
(3.2)
Здесь матрица A имеет блочно-диагональную структуру
A1
0
A
0
0
0
,
... A
0 ...
A 2 ...
...
0
причем ее блоки A1, …, A+ определяются следующим образом:
Ai = pi, если pi – вещественный полюс;
0
A i
pi pi 1
1
, если pi и pi+1 – комплексно сопряженные полюса.
pi pi 1
Матрица B представляет собой вектор-столбец, состоящий из нулей и
единиц, причем Bi = 0, если полюса pi и pi+1 комплексно сопряженные, в противном случае Bi = 1; матрица C есть вектор-строка, элементы которой также
определяются в зависимости от наличия мнимой части у соответствующего полюса передаточной функции:
Ci = ri, если pi – вещественный полюс;
Ci = –ripi+1 – ri+1pi, Ci+1 = ri + ri+1, если pi и pi+1 – комплексно сопряженные полюса.
Последовательно изменяя индекс i от 1 до n и анализируя полюса pi передаточной функции W(s), определим ненулевые элементы матриц A, B, C:
49
если pi – вещественный полюс:
Aii = pi; Bi = 1; Ci = ri;
если pi и pi+1 – комплексно сопряженные полюса:
Ai i+1 = 1, Ai+1 i = – pipi+1, Ai+1 i+1 = pi + pi+1;
Bi+1 = 1; Ci = –ripi+1 – ri+1pi, Ci+1 = ri + ri+1.
Систему уравнений (3.2) можно переписать несколько в ином виде, удобном для синтеза оптимального регулятора с использованием пакета Matlab:
x Ax Bu ξ,
z Cx υ ,
(3.3)
где x – многомерный белый шум с матрицей интенсивностей Nx = BBTN1.
Сформируем управление u(t) как выход линейно-квадратичного гауссовского регулятора (LQG-регулятора) в соответствии с критерием качества
J M (x T ()Qx () u 2 ()) d min,
0
(3.4)
где Q – задаваемая матрица весовых коэффициентов.
Оптимальный LQG-регулятор представляет собой комбинацию фильтра
Калмана и LQ-регулятора, синтезированного в соответствии с критерием (3.4)
[5, 15] и описываемого уравнением
u BHx̂,
где H – положительно определенное решение алгебраического уравнения Риккати [4]
A T H HA HBB T H Q 0.
50
В пакете Matlab процедура синтеза LQG-регулятора автоматизирована и
осуществляется с помощью встроенной функции lqg, аргументами которой являются матрицы, определяющие модель системы в пространстве состояний,
интенсивности входных шумов и заданные веса в критерии качества.
Матрица Q в квадратичном критерии, как правило, задается диагональной, что позволяет придать весовым коэффициентам Qjj достаточно ясный физический смысл, поскольку в этом случае
x T Qx Q11 x12 Q22 x22 ... Qnn xn2 .
А именно, увеличение того или иного весового коэффициента Qjj приведет к дополнительным ограничениям на возможность увеличения переменной
состояния xi по сравнению с другими переменными состояния и сигналом
управления. Тем не менее, в том случае, когда математическая модель объекта
управления получена без опоры на физический смысл переменных состояния,
выбор весовых коэффициентов может вызвать существенные трудности. Избежать этих трудностей позволит полученное выше описание упругого объекта, в
котором каждая из переменных состояния формирует одну из составляющих
движения объекта.
Пусть pi и pi+1 – комплексно сопряженные полюса передаточной функции
W(s). Тогда переменная состояния xi в первом приближении будет описывать
гармонические колебания с частотой i = Im pi и некоторой начальной фазой i:
xi (t ) xi max sin( i t i ),
где xi max – амплитуда колебаний. Соответственно, для переменной xi+1 max запишем
51
xi 1 (t ) xi max i cos( i t i ) xi 1max cos( i t i ).
Если колебания на частоте i являются преобладающими, выходной сигнал z
будет иметь вид
~
z (t ) Ci xi (t ) Ci 1 xi 1 (t ) (Ci sin( i t i ) Ci 1i cos( i t i )) xi max
C
i sin( i t i ) Ci 1 cos( i t i ) xi 1max ,
i
а его максимальное значение составит
Ci2
~
z max xi 1max
Ci21 .
2
i
Заметим, что передаточная функция от управления u(t) к переменной состояния xi(t) соответствует колебательному звену. Поэтому управление по переменной xi приведет к подавлению не только упругих колебаний на частоте i,
но и низкочастотной составляющей спектра входного сигнала, что в свою очередь вызовет увеличение перерегулирования в системе. Руководствуясь данными соображениями, зададим нулевыми весовые коэффициенты Qii при переменных xi (а также при переменных, формирующих апериодические составляющие движения) и ориентировочно определим коэффициенты Qi+1 i+1 при производных xi+1 = ẋi, воспользовавшись методом равных вкладов [17]. Согласно
данному методу весовые коэффициенты в критерии (3.4) следует выбирать таким образом, чтобы максимально допустимые отклонения переменных состояния вносили в функционал качества одинаковый вклад, равный вкладу максимально допустимого сигнала управления. Положив приблизительно равным
единице коэффициент передачи регулятора на частоте резонанса, получим
52
2
2
~
u max
z max
Ci2
Ci2
2
Qi 1 i 1 2
2
2 Ci 1
Ci21.
2
xi 1max xi 1max i
(Im pi )
(3.5)
Следует учесть, что значения коэффициентов, вычисленных по формуле
(3.5) являются ориентировочными и в ходе выполнения работы должны быть
скорректированы таким образом, чтобы обеспечить выполнение требований к
замкнутой системе по быстродействию (длительность переходного процесса не
должна превышать заданного значения tп) и запасу устойчивости (перерегулирование не должно превышать заданного значения ). При этом сигнал управления по абсолютной величине не должен превышать заданного значения umax.
3.2.Описание программы моделирования
Лабораторная работа выполняется в среде Matlab 6.5 с использованием
пакета Simulink.
Расчет регулятора и оценка качества переходного процесса выполняются
с использованием m-файла. Передаточная функция объекта определяется посредством задания ее числителя и знаменателя в виде символьных выражений.
После запуска файла на выполнение производится построение переходной и
амплитудно-частотной характеристик замкнутой системы.
Для исследования работы системы в условиях возмущений формируется
Simulink-модель в соответствии с рис. 3.3.
53
Band-Limited
White Noise
num(s)
1
den(s)
Constant
Transfer Fcn
Scope
x' = Ax+Bu
y = Cx+Du
Scope1
State-Space
Band-Limited
White Noise1
Рис. 3.3. Модель замкнутой системы в пакете Simulink
Возмущающие входные воздействия моделируются широкополосными
процессами с постоянной времени, меньшей периода колебаний наиболее высокочастотной колебательной составляющей движения объекта. Для всех вариантов интенсивности белых шумов принимаются равными N1 = 10–2 c,
N2 = 10–4 c.
Выходной сигнал системы можно пронаблюдать в блоке Scope, сигнал
управления – в блоке Scope1.
3.3.Порядок выполнения лабораторной работы
1. Получить вариант задания у преподавателя.
2. Запустить Matlab, создать новый mdl-файл.
3. Сформировать Simulink-модель в соответствии с рис. 3.3. Ввести:
в блок Constant
Constant value:
1
в блок Band-Limited White Noise
Noise power: N1
Sample time: 1E-4
в блок Transfer Fcn
Numerator:
num
54
Denominator: den
в блок State-Space
A:
af
B:
bf
C:
cf
D:
df
Initial conditions:
0
в блок Band-Limited White Noise1
Noise power: N2
Sample time: 1E-4
в блоки Scope, Scope1 (меню ′Scope′ parameters)
вкладка General:
Number of axes: 1
Time range:
auto
вкладка Data history:
Limit data points to last:
4. В меню Simulation выбрать пункт Simulation Parameters и установить
время моделирования 50 с.
5. Создать m-файл в соответствии с приложением 3. Ввести в соответствующие строки m-файла выражения для числителя и знаменателя передаточной функции в соответствии с вариантом задания (табл. 3.1).
55
Табл. 3.1 Варианты заданий
Вари, %
ант
1
60
2
50
3
30
4
0
5
60
6
50
7
60
tп, с
umax
50
10
5
20
25
50
20
2
2
1
2
3
2
4
W(s)
1
0,3( s 0,05)(s 250)( s 2 180s 3 10 4 )(s 2 90 s 2 10 4 )(s 2 150 s 5 10 4 )
( s 0,1)(s 2 10 s 2 10 5 )(s 2 5s 8 10 4 )(s 2 0,5 s 10 4 )(s 2 0,05 s 4)
2
0,15( s 0,5)( s 1)( s 100)( s 250)( s 2 500 s 10 5 )(s 2 30 s 4 10 4 )
( s 1,7)( s 2 10 s 2 10 5 )(s 2 5s 8 10 4 )(s 2 0,1s 600)( s 2 1,2 s 4)
3
0,3( s 1)( s 250)( s 2 500 s 10 5 )(s 2 0,5 s 100)( s 2 30 s 4 10 4 )
( s 1,7)( s 2 10 s 2 10 5 )(s 2 5s 8 10 4 )(s 2 0,1s 600)( s 2 1,2 s 4)
4
0,1( s 0,15)( s 250)( s 2 28 s 500)( s 2 85 s 2 10 4 )(s 2 300 s 4 10 4 )
( s 0,07)( s 2 3s 3 10 5 )(s 2 3s 10 4 )(s 2 0,5 s 100)( s 2 0,05 s 20)
5
0,2( s 0,1)( s 300)( s 2 70 s 1,5 10 4 )(s 2 150 s 4 10 4 )(s 2 400 s 5 10 4 )
( s 2)( s 2 7 s 10 5 )(s 2 1,5 s 5 10 4 )(s 2 0,6 s 5000)( s 2 0,1s 6)
6
0,2( s 0,05)( s 200)( s 2 140 s 2,2 10 4 )(s 2 75 s 8000)( s 2 30 s 4000)
( s 0,1)( s 2 s 10 5 )(s 2 2,5 s 9000)( s 2 0,5 s 850)( s 2 0,05 s 4)
7
0,15( s 0,1)( s 600)( s 2 50 s 1000)( s 2 350 s 4 10 4 )(s 2 900 s 5 10 4 )
( s 2)( s 2 5s 10 5 )(s 2 0,2 s 5 10 4 )(s 2 0,2 s 5000)( s 2 2s 6)
6. Нажатием клавиши [F5] запустить файл на выполнение.
7. Проанализировать построенные в раскрывшихся графических окнах
переходную (Step Response) и амплитудно-частотную (Bode Diagram) характе-
56
ристики. Выделить одну-две наиболее значимые моды колебаний, приблизительно определить их частоты.
8. Уточнить значения отмеченных в п. 7 частот колебаний, сопоставив их
с мнимыми частями значений полюсов pi передаточной функции W(s), отображаемых в командном окне Matlab. Зафиксировать индексы переменных состояния xi и xi+1, отвечающих за колебательные движения объекта с данными частотами.
9. В соответствии с методом равных вкладов ориентировочно найти по
формуле (3.5) значения весовых коэффициентов Qi+1 i+1 при переменных состояния xi+1, отмеченных в п. 8. Ввести полученные коэффициенты в m-файл.
Для нахождения мнимой части числа в Matlab используется функция
imag, например:
>> C(5)^2/imag(p(5))^2+C(6)^2
10. Запустить m-файл на выполнение. Раскрыть окно с построенными переходными характеристиками исходной (разомкнутой) системы (красная
штриховая линия) и скорректированной системы с регулятором в обратной
связи (синяя сплошная линия).
11. Щелчком правой кнопки мыши в области построения графика вызвать контекстное меню и выбрать Characteristics (характеристики), Peak Response (максимальное значение). Подвести курсор последовательно к выделенным точкам и определить перерегулирование (Overshoot) исходной и скорректированной систем.
12. В меню Characteristics выбрать пункт Setting Time (время переходного процесса). Подвести курсор последовательно к выделенным точкам и
определить длительность переходного процесса исходной и скорректированной систем.
13. Запустить mdl-файл на выполнение. Убедиться в отсутствии ярко выраженных колебаний на резонансных частотах в выходном сигнале z
(блок
Scope). Определить максимальное значение сигнала управления (блок Scope1).
57
14. Сравнить полученные в пп. 11-13 значения показателей качества системы управления с заданными в табл. 3.1.
15. В случае если синтезированная система по какому-либо показателю
не удовлетворяет заданным требованиям, изменить значения некоторых весовых коэффициентов (кроме последнего, равного 1) и повторить действия в соответствии с пп. 10-14. При корректировке весовых коэффициентов следует
руководствоваться следующими соображениями:
изменению подлежат, прежде всего, коэффициенты при производных
составляющих колебательного движения Qi+1 i+1;
пропорциональное увеличение весовых коэффициентов при переменных состояния приводит к улучшению качества переходного процесса, но увеличивает управление;
увеличение весового коэффициента при переменной состояния, формирующей одну из мод колебаний, по сравнению с весовым коэффициентом
при переменной состояния, формирующей другую моду, усиливает подавлению первой, но ослабляет подавление второй моды;
наибольшее влияние на показатели качества замкнутой системы, как
правило, оказывает коэффициент Qi+1 i+1 при переменной состояния xi+1, отвечающей за колебания системы с наиболее низкой из резонансных частот;
с целью уменьшения длительности переходного процесса можно задавать ненулевые весовые коэффициенты при переменных состояния, отвечающих за апериодическое движение, однако перерегулирование в этом случае
может сильно возрасти.
16. Определить передаточную функцию регулятора Wр(s), удовлетворяющего всем заданным требованиям. Для этого вызвать функцию tf(sys2) в командном окне Matlab. Записать полученное выражение.
58
3.4.Оформление отчета
Отчет о лабораторных исследованиях должен содержать:
структурную схему замкнутой системы;
исходные данные для синтеза;
значения весовых коэффициентов, использовавшиеся при оптимизации закона управления;
полученное выражение для передаточной функции регулятора;
графики переходной и амплитудно-частотной характеристик исходной
и скорректированной систем;
выводы по работе.
3.5.Контрольные вопросы
1. Как связаны полюса передаточной функции с частотами составляющих
колебательного движения?
2. Как изменится вид разложения передаточной функции на простейшие
дроби, если порядок ее числителя будет равен порядку знаменателя?
3. Все полюса передаточной функции замкнутой системы вещественные.
Обязательно ли переходная характеристика будет иметь апериодический характер? Почему?
4. Для чего в процессе синтеза закона управления передаточная функция
упругого объекта раскладывается на простейшие дроби?
5. Передаточная функция объекта с одним входом и одним выходом имеет m нулей и n полюсов. Каковы будут размерности матриц A, B, C при его
описании в пространстве состояний?
6. Что можно сказать о полюсах передаточной функции объекта, если его
матрица динамики диагональная?
7. Какова особенность использования квадратичного критерия оптимальности в стохастических системах?
59
8. С какой целью в работе используется метод равных вкладов? В чем его
суть?
9. Каким будет результат синтеза оптимального LQG-регулятора, если
все весовые коэффициенты при переменных состояния задать равными нулю?
10.Опишите способы формирования линейных стационарных объектов в
среде Matlab.
60
4. СИНТЕЗ АЛГОРИТМА ЭКСТРАПОЛЯЦИИ УЗКОПОЛОСНОГО
СЛУЧАЙНОГО СИГНАЛА
Цель работы: изучение процедуры экстраполяции случайных процессов
с применением алгоритма идентификации.
4.1.Методические указания по подготовке к работе
Управление движением морских судов происходит в условиях интенсивных волновых возмущений. Эффективным путем уменьшения влияния морского волнения на качество управления является использование как можно более
полной информации об угловом движении судна. В частности, возможность
получения экстраполированных (упрежденных) значений углов качки может
способствовать предотвращению аварийной ситуации в работе движителя судна, когда ходовой винт при определенных углах качки и уровне волнения может оказаться на некоторое время вне водной среды. Кроме того, решение задачи экстраполяции (прогнозирования) качки требуется при работе систем
управления посадкой корабельной авиации, когда возникает необходимость
определения взаимного расположения палубы корабля и летательного аппарата.
Возможные подходы к решению задачи прогнозирования могут быть основаны на построении оптимальных экстраполяторов с применением теории
винеровской или калмановской фильтрации [8, 18, 19]. Общий недостаток получаемых решений связан с априорным представлением качки в виде стационарного случайного процесса с известными характеристиками. В то же время
известно, что параметры качки (дисперсия, преобладающая частота, ширина
спектра) зависят от целого ряда факторов, таких как балльность морского волнения, величина курсового угла к волне, скорость движения судна, направление ветра и т.д. [20]. Отсутствие точных данных о параметрах качки может
привести к значительным ошибкам в прогнозировании. Для получения досто61
верных оценок параметров качки с высоким быстродействием может быть использована идентификационная процедура, основанная на нелинейных алгоритмах обработки измерений. Синтез таких алгоритмов удобно проводить в
предположении о возможности представления процесса качки судна на ограниченных отрезках времени (порядка нескольких периодов) в виде гармонических колебаний, происходящих с частотой, близкой к преобладающей [21]. При
таком описании предполагается, что отклонения амплитуды и частоты качки от
средних значений являются медленноменяющимися случайными процессами с
корреляционными функциями, близкими к экспоненциальным.
Считая несущественными изменения амплитуды колебаний по сравнению с изменениями частоты, запишем
(t ) A0 sin (0 (t ))t 0 ,
(4.1)
где (t) – текущее значение угла качки, A0 – амплитуда качки, 0 – среднее значение частоты качки, (t) – случайные девиации частоты, 0 – случайная
начальная фаза.
Дважды продифференцировав выражение (4.1) по времени с учетом малости производной ̇, получим систему дифференциальных уравнений относительно переменных состояния x1 = и x 2 :
x1 x2 ,
(4.2)
x 2 (0 x3 ) 2 x1 ,
(4.3)
где x3 = . Дополним полученную систему уравнением формирующего фильтра для экспоненциально коррелированного процесса (t) [14]:
x3 x3 2x,
(4.4)
62
где – показатель затухания корреляционной функции процесса (t),
– среднеквадратическое отклонение частоты качки от среднего значения,
x(t) – входной белый шум формирующего фильтра.
В векторно-матричной форме уравнения (4.2)-(4.4) будут иметь вид
x f (x) Bx,
T
T
где x – вектор состояния, B 0 0 2 – матрица воз-
мущений.
Компоненты вектор-функции f(x) составят
f1 (x) x2 ,
(4.5)
f 2 (x) (0 x3 ) 2 x1 ,
(4.6)
f 3 (x) x3 2x.
(4.7)
Предположим, что измерению подвергаются угол качки и скорость его
изменения. При учете случайных (белошумных) погрешностей измерителей
уравнения измерений представляются в виде
z1 x1 1 ,
(4.8)
z 2 x2 2 ,
(4.9)
или в векторно-матричной форме
z Hx v,
63
где z = (z1 z2)T – вектор измерений, v = (
1
2)
T
– вектор шумов измерений,
1 0 0
– матрица наблюдений.
H
0 1 0
Поскольку функция f(x) является нелинейной, для решения задачи оценивания компонент вектора состояния будем использовать нелинейный фильтр
Калмана-Бьюси в виде [14]
xˆ f (xˆ ) K (z Hxˆ ),
(4.10)
f (xˆ ) 1
R ,
x T
(4.11)
KP
f (xˆ )
f (xˆ )
P
PP
PH T R 1HP BB T ,
T
x
x
T
(4.12)
где x̂ – текущая оценка вектора состояния; P – матрица ковариаций вектора
ошибок ε x x̂; R – диагональная матрица интенсивностей шумов измерений
(с заданными компонентами R1, R2); K – матричный коэффициент усиления
фильтра.
В соответствии с выражениями (4.5)-(4.7) производная
f (xˆ )
будет иметь
x T
вид
f1
x1
f (x) f 2
x T x1
f 3
x
1
f1
x2
f 2
x2
f 3
x2
f1
x3
0
f 2
(0 x3 ) 2
x3
0
f 3
x3
0 2 x1 (0 x3 ) .
0
1
0
(4.13)
Начальными условиями для интегрирования уравнений (4.10)-(4.12) являются априорная оценка вектора состояния xˆ (0) 0 и априорная матрица ко64
вариаций вектора ошибок P(0) = P0, представляющая собой диагональную матрицу, содержащую средние квадраты компонент вектора состояния:
A02 2
P0 0
0
0
0
A02 (02 2 ) 2 0 .
0
2
Для получения значения угла качки, упрежденного на время , введем в
рассмотрение линеаризованную систему уравнений (4.2), (4.3) относительно
оценок переменных состояния xˆ1 , xˆ 2 :
xˆ1 xˆ 2 ,
(4.14)
ˆ ) 2 xˆ1 ,
xˆ 2 (0
(4.15)
ˆ xˆ3 – текущая оценка девиации частоты, и воспользуемся формулой
где
Коши при нулевом входном воздействии:
xˆ (t ) Φ(t , t )xˆ (t ),
0
где (t + , t) = eA – фундаментальная матрица, A
ˆ )2
(0
1
– мат0
рица динамики системы (4.14), (4.15). Матричная экспонента может быть
найдена с помощью преобразования Лапласа [3]:
e A L1 ( pE A) 1 .
Для рассматриваемой системы получим
65
p
pE A
ˆ )2
(0
1
;
p
p
2
ˆ )2
p (0
1
( pE A)
(0
ˆ )2
2
ˆ )2
p (0
e
A
1
ˆ )2
p 2 (0
;
p
ˆ )2
p 2 (0
1
ˆ )
ˆ )
cos( 0
sin( 0
ˆ
0
.
(
ˆ ) sin( 0
ˆ )
ˆ )
cos( 0
0
Отсюда выражение для экстраполированной оценки угла качки примет
вид
ˆ (t ) xˆ1 (t ) xˆ1 (t ) cos( 0 xˆ3 (t ))
xˆ 2 (t )
sin( 0 xˆ3 (t )).
0 xˆ3 (t )
При этом текущие оценки компонент вектора состояния xˆ j (t ) определяются путем интегрирования матричных уравнений (4.10)-(4.12).
4.2.Описание программы моделирования
Лабораторная работа выполняется в среде Matlab 6.5 с использованием
пакета Simulink.
Ввод исходных данных и расчет матриц B, R, P0 производятся с использованием m-файла.
Для исследования работы системы в условиях возмущений формируется
mdl-файл в соответствии с рис. 4.1. Модель фильтра Калмана-Бьюси (блок
Kalman Filter) имеет вид, приведенный на рис. 2.4, при этом его внутренние
блоки MATLAB Function характеризуются функциями (4.5)-(4.9), (4.13). Вычисление
упрежденных
значений
углов
качки
осуществляется
в
блоке
Extrapolation.
66
При описании полезного сигнала используется модель в виде гармонических колебаний с заданной частотой 1 и начальной фазой 1:
(t ) A0 sin 1t 1 .
При этом предполагается, что частота 1 может существенно отличаться от частоты настройки фильтра 0.
Возмущающие входные воздействия моделируются широкополосными
процессами с малой постоянной времени (0,01 с).
В блоке Scope1 можно пронаблюдать полезный сигнал, его оценку, приведенную ко времени экстраполяции, а также ошибку оценивания. В блоке
Scope формируются графики изменения во времени элементов ковариационной
матрицы P.
Band-Limited
White Noise
u
z
Sine Wave
x^
MATLAB
Function
P
Extrapolation
Kalman Filter
Scope
Sine Wave1
Band-Limited
White Noise1
theta, rad
R2D
theta, deg
Radians
to Degrees
Transport
Delay
theta_error, rad
R2D
theta_error, deg
Radians
to Degrees1
Scope1
Рис. 4.1. Модель исследуемой системы в пакете Simulink
67
4.3.Порядок выполнения лабораторной работы
1. Получить вариант задания у преподавателя.
2. Запустить Matlab, создать новый mdl-файл.
3. Сформировать Simulink-модель в соответствии с рис. 4.1. Ввести:
в блок Sine Wave
Amplitude:
Bias:
A0
0
Frequency (rad/sec): omeg1
Phase (rad):
fi1
Sample time: 0
в блок Sine Wave1
Amplitude:
Bias:
A0*omeg1
0
Frequency (rad/sec): omeg1
Phase (rad):
fi1+pi/2
Sample time: 0
в блок Band-Limited White Noise
Noise power: R1
Sample time: 1E-2
в блок Band-Limited White Noise1
Noise power: R2
Sample time: 1E-2
в блок Extrapolation
MATLAB Function:
u(1)*cos((omeg0+u(3))*tau)+u(2)/(omeg0+u(3))*sin((omeg0+u(3))*tau)
в блок Transport Delay
Time delay:
tau
в блок Scope1 (меню ′Scope′ parameters)
вкладка General:
68
Number of axes: 2
Time range:
20
вкладка Data history:
Limit data points to last:
4. Раскрыть блок Kalman Filter. Сформировать Simulink-модель в соответствии с рис. 2.4. Ввести:
в блок Constant
Constant value:
R^-1
в блоки Product, Product1, Product2, Product3
Multiplication:
Matrix(*)
в блок Integrator
Initial condition:
[0; 0; 0]
в блок f(x^,u)
MATLAB Function: [u(3); -(omeg0+u(4))^2*u(2); -alf*u(4)]
Output dimensions: 3
в блок h(x^)
MATLAB Function: [u(1); u(2)]
Output dimensions: 2
в блок df(x^,u)/dx'
MATLAB Function:
[0; -(omeg0+u(4))^2; 0; 1; 0; 0; 0; -2*u(2)*(omeg0+u(4)); -alf]
Output dimensions: 9
в блок Integrator1
Initial condition:
P0
в блок dh(x^)/dx'
MATLAB Function: [1; 0; 0; 1; 0; 0]
Output dimensions: 6
в блок Reshape
Output dimensionality:
Customize
69
Output dimensions: [3,3]
в блок Reshape1
Output dimensionality:
Customize
Output dimensions: [2,3]
в блоки Math Function, Math Function1
Function:
transpose
в блок Constant1
Constant value:
B*B'
5. В меню Simulation выбрать пункт Simulation Parameters и установить
время моделирования 20 с.
6. Создать m-файл в соответствии с приложением 4. Ввести значения параметров в соответствии с вариантом задания (табл. 4.1), обращая внимание на
их размерности. Задать 1 = 0,5 с–1, = 1 с.
Табл. 4.1. Варианты заданий
ВариA0, град 1, град
ант
1
3
0
2
4
45
3
5
60
4
5
90
5
6
225
6
7
180
7
8
270
R1,
град2c
10–6
10–5
10–6
10–5
10–6
10–5
10–6
R2,
град2/c
10–4
210–4
510–5
210–4
10–4
210–4
510–5
0, с–1
, с–1
, с–1
3
3
3
2
3
3
2
0,2
0,2
0,3
0,3
0,3
0,2
0,2
1/300
1/360
1/200
1/100
1/180
1/150
1/240
7. Нажатием клавиши [F5] запустить m-файл на выполнение.
8. Вернуться в окно созданной Simulink-модели. Нажатием кнопки Start
Simulation запустить mdl-файл на выполнение.
9. Двойным щелчком мыши по блоку Scope1 раскрыть окно со строящимися графиками:
полезного сигнала и его экстраполированной оценки ̂;
ошибки экстраполяции ̂.
70
10.Определить максимальную ошибку экстраполяции m в установившемся режиме.
11.Повторить процедуру моделирования, задавая значения 1 и в соответствии с табл. 4.2. Зафиксировать полученные значения ошибок экстраполяции. Проанализировать зависимость точности экстраполяции от частоты полезного сигнала.
Табл. 4.2. Максимальная ошибка экстраполяции m, град
0,5
1
2
3
4
1, с–1
1
, с 3
5
5
6
4.4.Оформление отчета
Отчет о лабораторных исследованиях должен содержать:
математическую модель динамики системы;
модель реального сигнала;
заполненную табл. 4.2;
выводы по работе.
4.5.Контрольные вопросы
1. Каковы достоинства и недостатки алгоритмов прогнозирования, основанных на теории оптимальной фильтрации?
2. Каким образом может быть описана качка судна на нерегулярном волнении?
3. Какие параметры характеризуют качку судна?
4. Для чего в рассмотренном алгоритме экстраполяции используется
идентификационная процедура?
5. Почему действительная матрица ковариаций вектора ошибок может
отличаться от расчетной, вырабатываемой в фильтре Калмана?
71
6. Какой вид имеет общее решение уравнения x 02 x 0 ? Какие начальные условия необходимо задать, чтобы найти его частное решение?
7. Что такое фундаментальная матрица линейной системы? При решении
каких задач она используется?
8. Как найти фундаментальную матрицу линейной стационарной системы?
9. Укажите причины возрастания ошибок экстраполяции при увеличении
времени упреждения.
10. Какое соотношение амплитуд и какую разность фаз имеют гармонический сигнал и его производная?
72
ЛИТЕРАТУРА
1. Андриевский Б.Р., Фрадков А.Л. Избранные главы теории автоматического управления с примерами на языке MATLAB. – СПб.: Наука, 1999. –
467 с.
2. Андриевский Б.Р., Фрадков А.Л. Элементы математического моделирования в программных средах MATLAB 5 и Scilab. СПб.: Наука, 2001. – 50 с.
3. Методы классической и современной теории автоматического управления: Учеб. в 3-х т. Т.1: Анализ и статистическая динамика систем автоматического управления / Под ред. Н.Д. Егупова. – М.: Изд-во МГТУ им. Н.Э. Баумана, 2000. – 748 с.
4. Методы классической и современной теории автоматического управления: Учеб. в 3-х т. Т.2: Синтез регуляторов и теория оптимизации систем автоматического управления / Под ред. Н.Д. Егупова. – М.: Изд-во МГТУ им.
Н.Э. Баумана, 2000. – 736 с.
5. Методы классической и современной теории автоматического управления: Учеб. в 3-х т. Т.3: Методы современной теории автоматического управления / Под ред. Н.Д. Егупова. – М.: Изд-во МГТУ им. Н.Э. Баумана, 2000. –
748 с.
6. Сергиенко А.Б. Цифровая обработка сигналов. СПб.: Питер, 2003. –
604 с.
7. Ерофеев А.А. Теория автоматического управления: Учеб. для вузов. –
СПб.: Политехника, 2002. – 302 с.
8. Радиоавтоматика: Учеб. пособие / Под ред. В.А. Бесекерского. – М:
Высш. школа, 1985. – 271 с.
9. Вагущенко Л.Л., Цымбал Н.Н. Системы автоматического управления
движением судна. – Одесса: ЛАТСТАР, 2002. – 310 с.
10.Березин С.Я., Тетюев Б.А. Системы автоматического управления движением судна по курсу. – Л.: Судостроение, 1990. – 256 с.
73
11.Соболев Г.В. Управляемость корабля и автоматизация судовождения.
– Л.: Судостроение, 1976. – 480 с.
12.Анучин О.Н., Емельянцев Г.И. Интегрированные системы ориентации
и навигации для морских подвижных объектов. – СПб.: ГНЦ РФ – ЦНИИ
"Электроприбор", 2003. – 390 с.
13.Бромберг П.В. Теория инерциальных систем навигации. – М.: Наука,
1979. – 278 с.
14.Степанов О.А. Применение теории нелинейной фильтрации в задачах
обработки навигационной информации. – СПб: ГНЦ РФ – ЦНИИ "Электроприбор", 1998. – 370 с.
15.Брайсон А., Хо Ю-Ши. Прикладная теория оптимального управления:
Пер. с англ. под ред. А.М. Летова. – М.: Мир, 1972. – 544 с.
16.Бродский С.А., Небылов А.В., Панферов А.И. Контроль параметров и
стабилизация упругих колебаний сложных объектов и конструкций // Гироскопия и навигация. – 2006. – №4(55). – С. 105
17.Колганов А.Р. Основные разделы современной теории автоматического
управления:
Электронный
конспект
лекций.
–
http://elib.ispu.ru/library/lessons/kolganov2/index.html
18.Ривкин С.С. Прогнозирование положения палубы корабля на качке //
Вопросы кораблестроения, Сер. Навигация и гироскопия. – 1980. – Вып. 49. –
С. 55-64.
19.Крайнов В.И., Тупысев В.А. Об экстраполяции вырабатываемых значений углов качки корабля // Гироскопия и навигация. – 1994. – №1(4). – С. 5864.
20.Бородай И.К., Нецветаев Ю.А. Качка судов на морском волнении. – Л.:
Судостроение, 1969. – 432 с.
21.Максимова Г.Н. Экстраполирование узкополосного сигнала // Изв. вузов. Приборостроение. – 1973.
74
ПРИЛОЖЕНИЕ 1
Текст m-файла к лабораторной работе №1
%Начальные приближения для коэффициентов ПИД-регулятора
kpid = [-1
-0.01
-1];
kp = kpid(1);
ki = kpid(2);
kd = kpid(3);
%Параметры модели судна
k1 = -0.13; % 1/c
tau1 = -60; % c
tau2 = 6;
% c
tau3 = 15;
% c
c2 = 0;
% c
c3 = -700;
% c^2
b1 = k1*tau3/tau1/tau2;
b2 = k1*(tau1*tau2 -tau1*tau3 - tau2*tau3)/tau1^2/tau2^2;
%Параметры рулевого привода
dmin = 0.4*pi/180; %граница зоны нечувствительности
%привода, рад
dmax = 35*pi/180; %максимальный угол перекладки руля, рад
ddmax = 4*pi/180; %максимальная скорость перекладки руля,
%рад/с
tau4 = 0.3;
%постоянная времени привода, с
K3 = 90*pi/180; %отклонение курса заданного, рад
75
ПРИЛОЖЕНИЕ 2
Текст m-файла к лабораторной работе №2
clear; %очистка рабочей области
close all; %закрытие окон
rand('state',0); %установка датчика случайных чисел в
%исходное состояние
g = 9.81; %ускорение силы тяжести
%СКО угловой скорости дрейфа, 1/с
sig1 = 0.005/180*pi;
sig2 = 0.005/180*pi;
sig3 = 0.005/180*pi;
%постоянные времени угловой скорости дрейфа, с
T1 = 10;
T2 = 10;
T3 = 10;
%СКО составляющих абсолютного ускорения, м/с^2
siggx = 0.005;
siggy = 0.005;
siggz = 0.005;
%постоянные времени составляющих абсолютного ускорения, c
Tgx = 1;
Tgy = 1;
Tgz = 1;
B(7,1) = sig1*sqrt(2/T1);
76
B(8,2) = sig2*sqrt(2/T2);
B(9,3) = sig3*sqrt(2/T3);
B(10,4) = siggx*sqrt(2/Tgx);
B(11,5) = siggy*sqrt(2/Tgy);
B(12,6) = siggz*sqrt(2/Tgz);
%амплитуды составляющих угловой скорости, 1/c
Omeg_x = 10*pi/180;
Omeg_y = 15*pi/180;
Omeg_z = 20*pi/180;
%частоты изменения составляющих угловой скорости, 1/с
Omegx = 2*pi/10;
Omegy = 2*pi/10;
Omegz = 2*pi/10;
%начальная фаза изменения угловой скорости
fi = unifrnd(-pi,pi);
%матрица интенсивностей шумов измерений
R = diag([10^-7; 10^-7; 10^-7]);
%начальные значения углов ориентации
tet0 = -1*pi/180;
psi0 = 1*pi/180;
gam0 = -1*pi/180;
x0 = [sin(tet0); cos(tet0)*cos(gam0);
-cos(tet0)*sin(gam0); -sin(psi0)*cos(tet0);
cos(psi0)*sin(gam0)+sin(psi0)*sin(tet0)*cos(gam0);
cos(psi0)*cos(gam0)-sin(psi0)*sin(tet0)*sin(gam0);
normrnd(0,sig1); normrnd(0,sig2); normrnd(0,sig3);
77
normrnd(0,siggx); normrnd(0,siggy); normrnd(0,siggz)];
%начальная оценка вектора состояния
xo0 = [0; 1; 0; 0; 0; 1; 0; 0; 0; 0; 0; 0];
P0 = diag([10^-2; 10^-2; 10^-2; 10^-2; 10^-2; 10^-3;
sig1; sig2; sig3; siggx; siggy; siggz].^2);
78
ПРИЛОЖЕНИЕ 3
Текст m-файла к лабораторной работе №3
clear; %очистка рабочей области
close all; %закрытие окон
N1 = 1E-2; %интенсивность входного шума
N2 = 1E-4; %интенсивность шума измерений
%определение коэффициентов числителя и знаменателя
%передаточной функции
syms 's';
num=sym2poly(expand(0.15*(s+0.1)*(s+600)*(s^2+50*s+1000)*
(s^2+350*s+40000)*(s^2+900*s+50000))); %числитель
den=sym2poly(expand((s+2)*(s^2+5*s+1E+5)*(s^2+0.2*s+5E+4)
*(s^2+0.2*s+5000)*(s^2+2*s+6))); %знаменатель
%разложение передаточной функции на простейшие дроби
[r,p,k] = residue(num,den)
sys1 = ss(tf(num,den));
%формирование матриц A, B, C
C = zeros(size(den)-[0 1]);
B = C';
for i = 1:size(p)
if imag(p(i))==0 %если полюс вещественный
A(i,i) = p(i);
C(i) = r(i);
B(i) = 1;
elseif B(i)==0 %если полюс комплексный, а предыдущий
%был вещественным
79
A(i,i+1) = 1;
A(i+1,i) = -p(i)*p(i+1);
A(i+1,i+1) = p(i)+p(i+1);
C(i) = -r(i)*p(i+1)-r(i+1)*p(i);
C(i+1) = r(i)+r(i+1);
B(i+1) = 1;
end
end
%из-за ошибок округления могут появиться мнимые части,
%следует их отбросить
A = real(A);
C = real(C);
%матрица интенсивностей белых шумов
V = [B*B'*N1 zeros(size(B),1); zeros(1,size(B)) N2];
%матрица весовых коэффициентов в критерии оптимальности
W = diag([0 0 0 0 0 0 0 0 0 1]);
[af,bf,cf,df] = lqg(A,B,C,0,W,V); %расчет LQG регулятора
sys2 = ss(af,bf,cf,df);
sys = sys1/(1+sys1*sys2); %замыкание обратной связи
%Построение графиков
%переходная характеристика разомкнутой системы
step(sys1,'r--');
hold on; %в том же окне
%переходная характеристика замкнутой системы
step(sys);
figure; %новое окно
80
bodemag(sys1,'r--'); %ЛАХ разомкнутой системы
hold on; %в том же окне
bodemag(sys); %ЛАХ замкнутой системы
81
ПРИЛОЖЕНИЕ 4
Текст m-файла к лабораторной работе №4
clear;
tau = 1; %время прогнозирования, c
omeg1 = 0.5; %частота качки, 1/c
A0 = 8*pi/180; %амплитуда качки, рад
fi1 = 3*pi/2; %начальная фаза, рад
R1 = 1E-6*(pi/180)^2; %интенсивность шума измерений
%угла качки, рад^2*с
R2 = 5E-5*(pi/180)^2; %интенсивность шума измерений
%угловой скорости, рад^2/c
omeg0 = 2; %центральная частота настройки фильтра, 1/с
sigomeg = 0.2; %СКО частоты, 1/c
alf = 1/240; %показатель затухания корреляционной функции
%смещения частоты, 1/c
R = [R1 0; 0 R2]; %матрица интенсивностей шумов измерений
B(3,1) = sigomeg*sqrt(2*alf); %матрица возмущений
P0 = diag([A0^2/2; A0^2*(omeg0^2+sigomeg^2)/2;
sigomeg^2]); %априорная матрица ковариаций
82
