close

Вход

Забыли?

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

?

Panferov 05B50E71DB

код для вставкиСкачать
Федеральное агенТство по образованию
Государственное образовательное учреждение
высшего профессионального образования
Санкт-петербургский государственный университет
аэрокосмического приборостроения
А. И. Панферов, А. В. Лопарев
КОМПЬЮТЕРНЫЙ АНАЛИЗ И
СИНТЕЗ СИСТЕМ ОРИЕНТАЦИИ,
СТАБИЛИЗАЦИИ И НАВИГАЦИИ
Учебно-методическое пособие
Санкт-Петербург
2008
УДК 004.9
ББК 32.97
П16
Рецензенты:
кафедра инерциальных навигационных систем Санкт-Петербургского
государственного университета информационных технологий, механики
и оптики; кандидат технических наук А. Е. Елисеенков
(ГНЦ РФ – ЦНИИ «Электроприбор»)
Утверждено
редакционно-издательским советом университета
в качестве учебного пособия
Панферов А. И., Лопарев А. В.
П00 Компьютерный анализ и синтез систем ориентации, стабилизации и навигации: учебное пособие / А. И. Панферов,
А. В. Лопарев. – СПб.: ГУАП, 2008. – 64 с.: ил.
ISBN 978-5-8088-0409-0
Учебно-методическое пособие содержит описание лабораторных
работ по исследованию динамических и точностных свойств систем
ориентации, стабилизации и навигации. Каждая работа содержит
методические указания к выполнению заданий, описание программы моделирования, контрольные вопросы. Приводятся необходимые
сведения из теории фильтрации и теории управления. В качестве среды моделирования используется пакет Matlab 6.5.
Пособие предназначено для студентов технических вузов, обучающихся по специальности 181200.
УДК 004.9
ББК 32.97
ISBN 978-5-8088-0409-8
© ГУАП, 2008
© А. И. Панферов,
А. В. Лопарев,
2008
ПРЕДИСЛОВИЕ
Развитие средств вычислительной техники предоставляет все
больше возможностей при решении задач теории автоматического
управления. Вычислительные устройства используются как на этапе синтеза при предварительном расчете характеристик системы,
так и непосредственно в контуре управления при ее практической
реализации.
Одним из наиболее популярных инструментов научных исследований и математических расчетов стала система Matlab, построенная на расширенном представлении и применении матричных
операций. Матрицы широко применяются при построении и математическом моделировании динамических систем и объектов, они
являются основой составления и решения уравнений состояния.
Однако в настоящее время Matlab далеко вышла за пределы специализированной матричной системы и является одним из наиболее
мощных универсальных средств компьютерных расчетов и моделирования. При этом большое внимание уделяется наглядности и
удобству в использовании системы, что наглядно проявляется в пакете расширения Matlab – Simulink.
Обладая широкими возможностями по созданию различных
демонстрационных примеров, система Matlab������������������
������������������������
нашла свое применение и в учебном процессе. Изданная учебно-методическая литература по анализу и синтезу динамических систем, так или иначе
опирающаяся на ����������������������������������������������
Matlab����������������������������������������
[1–6], охватывает широкий спектр теоретических и практических вопросов, однако исследовательскому аспекту в подобных работах отведена второстепенная роль. Вместе с
тем, подготовка специалистов и магистров по направлению 652300
«Системы управления движением и навигация» в высших учебных
заведениях подразумевает приобретение студентами практических
навыков проектирования сложных динамических систем. Такие
навыки могут быть получены в ходе выполнения лабораторных работ с применением средств вычислительной техники, и предлагаемое учебное пособие может восполнить недостаток учебного материала в данной области.
Тематика лабораторных работ охватывает вопросы построения
систем с заданным качеством переходного процесса, синтеза оптимальных и субоптимальных алгоритмов управления и обработки
информации по критерию наивысшей точности, анализа чувствительности систем при изменении параметров входных воздействий.
К каждой работе даны методические рекомендации по выполнению
исследований, подробное описание программы моделирования,
контрольные вопросы.
3
1. СИНТЕЗ СИСТЕМЫ СТАБИЛИЗАЦИИ СУДНА ПО КУРСУ
С ЗАДАННЫМ КАЧЕСТВОМ ПЕРЕХОДНОГО ПРОЦЕССА
Цель работы: изучение процедуры параметрического синтеза
законов управления с заданием допустимого «коридора» для переходной характеристики.
1.1. Методические указания по подготовке к работе
Как известно [1, 3, 7, 8], все реальные системы автоматики нелинейны. Физическим звеньям свойственны явления насыщения,
ограничения, люфта, гистерезиса и т.д. Теория нелинейных систем сложна, недостаточно хорошо развита, и, по существу, анализ
систем всегда является приближенным. В то же время существует
целый класс нелинейных систем, которые при определенных допущениях можно линеаризовать, т.е. рассматривать линейные математические модели вместо нелинейных. Для таких систем нелинейные свойства не играют существенной роли, и их исследование,
как правило, не представляет особых технических трудностей.
С другой стороны, существует обширное множество систем, для
которых нелинейные свойства являются принципиально важными и применение линейных моделей приводит к качественно неверным результатам. Причем сложность и большое разнообразие
свойств таких систем не позволяют использовать какой-либо универсальный подход к их анализу и синтезу. В связи с этим при исследовании нелинейных систем широкое распространение получили численные методы, которые могут применяться, в том числе, и
с целью их оптимизации.
Настоящая лабораторная работа посвящена синтезу регулятора
для нелинейной системы управления с заданным качеством переходного процесса. Данный подход к синтезу основывается на целенаправленном поиске таких значений настраиваемых параметров
в законе управления, чтобы траектории движения динамической
системы не выходили за пределы заданных областей («коридоров»)
[7]. Существо подхода заключается в минимизации некоторой целевой функции, определяющей меру выхода исследуемой траектории за пределы «коридора». Если определить эту функцию таким образом, чтобы ее нулевое значение соответствовало касанию
траектории границ допустимых областей, то разрешимость задачи
синтеза будет определяться наличием неположительных значений
функции.
Описанный подход к параметрическому синтезу законов управления реализован в математической среде Matlab Simulink в форме
4
пакета прикладных программ NCD-blockset. Это инструментальное
средство позволяет задавать границы допустимого «коридора» в
визуальном режиме, обеспечивает вычисление меры выхода за его
пределы, осуществляет запуск численного метода решения оптимизационной задачи и находит ее решение, если оно существует.
Рассмотрим методику использования пакета NCD-blockset на
примере синтеза системы стабилизации морского судна по курсу
(авторулевого). Укрупненная структурная схема системы представлена на рис. 1.1, на котором обозначены: УУ – устройство
управления (регулятор); ИУ – исполнительное устройство (рулевой
привод); ОУ – объект управления (судно); Kз – заданное значение
курса; K – текущее значение курса; ψ = Kз – K – отклонение курса
от заданного значения (ошибка регулирования); u – сигнал управления; β – угол перекладки руля.
Для упрощенного описания движения судна (без учета внешних
возмущений) будем использовать нелинейное дифференциальное
уравнение второго порядка в форме Номото [9, 10]:
 + (τ1 + τ2 )ω
 + ω + c2 ω ω + c3 ω3 = k1β + k1τ3β , τ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).
ψ
Kз
+
УУ
u
ИУ
β
ОУ
K
–
Рис. 1.1. Структурная схема авторулевого
5
Характер диаграммы поворотливости зависит от устойчивости судна на курсе, определяеω
мой знаками коэффициентов
в уравнении (1.1). Когда суд2
но теоретически устойчиво на
курсе, все коэффициенты в
1
этом уравнении положительω0
βп
ны (кривая 1). Если судно не
β
обладает устойчивостью прямолинейного движения, то
τ1 < 0, k1 < 0, c2 ≤ 0, а τ2 > 0,
τ3 > 0 (кривая 2). Знак c3  зависит не только от устойчивости
Рис. 1.2. Диаграммы поворотливости судна, но и от вида диаграммы
поворотливости при больших
перекладках руля [9]. Из опыта судовождения следует [11], что при
c2 ≤ 0 значения |ω | < |ω0| у судов неустойчивы и могут быть получены только в результате некоторого процесса регулирования движения. Поэтому на рис. 1.2 кривая на данном участке изображена
штриховой линией. При этом изменение угловой скорости на предельных углах обратной поворотливости ±βп происходит не плавно,
а скачком, в результате чего зависимость ω(β) приобретает вид петли гистерезиса.
Для анализа математической модели движения судна по курсу
перейдем от дифференциального уравнения (1.1) к описанию системы в пространстве состояний:
x1 = x2 + b1β,
x2 = −
c
c
τ1 + τ2
1
x2 −
x1 − 2 x1 x1 − 3 x13 + b2β, τ1τ2
τ1τ2
τ1τ2
τ1τ2
(1.3)
x3 = x1.
В данной модели угловой скорости ω соответствует переменная
состояния x1, а углу ψ – переменная состояния x3. Коэффициенты
b1,2 определяются из соотношений [1]:
b1 =
k1τ3
,
τ1τ2
b2 = k1
6
τ1τ2 − τ1τ3 − τ2 τ3
τ12 τ22
.
Работа системы стабилизации курса считается удовлетворительной, когда судно достаточно быстро и с малым перерегулированием
выходит на заданный курс и стабилизируется на нем с высокой точностью удержания (малым рысканием). Обеспечить выполнение
требований быстродействия и малой колебательности возможно
при использовании пропорционально-интегрально-дифференциального закона управления (ПИД-регулятора) [9–11]:
u = kп ψ + kд ψ + kи ∫ ψdt, (1.4)
где kп, kд, kи – коэффициенты регулятора. При этом управление по
интегралу от углового рассогласования обеспечивает нулевую статическую ошибку регулирования, а управление по производной
служит для демпфирования собственных колебаний системы [7,
9].
Выражению (1.4) соответствует передаточная функция регулятора
WУУ (s) = kп + kд s + kи s.
Поскольку использование интегральной составляющей в законе
управления не способствует улучшению качества переходного процесса, ее целесообразно включать только в режиме стабилизации.
Поэтому при отклонениях курса от заданного значения на величину 5° и более будем полагать kи = 0.
Управление u отрабатывается рулевым приводом, математическую модель которого примем в виде [11]
ε = u − β,
ε = (ε − δmin sign ε) ⋅ Φ( ε − δmin ), β = sat( ε τ4 , δ max ) ⋅ Φ(δmax − β ),
(1.5)
где δmin – минимальный угол перекладки руля, отрабатываемый
приводом (половина ширины зоны нечувствительности); δmax – максимальный угол перекладки руля; δ max – максимальная скорость
перекладки руля; τ5 – постоянная времени привода; sat(x, xmax) –
функция насыщения:
x, x ≤ xmax ,

sat(x, xmax ) = xmax , x > xmax ,
 −x , x < −x ;
max
 max
Φ(x) – единичная функция Хевисайда:
7
1, x ≥ 0,
Φ (x) = 
0, x < 0.
При δmin ≤ u ≤ δmax, u ≤ δ max система (1.5) допускает линеаризованное представление:
β = (u − β) τ ,
4
что позволяет рассматривать привод как апериодическое звено первого порядка с передаточной функцией
WРП (s) =
1
.
1 + τ 4s
Однако в том случае, когда указанные условия не выполняются, существенную роль играют присущие рулевому приводу нелинейности: насыщение и нечувствительность. Первое из указанных
свойств привода проявляется в переходном режиме системы управления, второе – в режиме стабилизации.
Зададим «коридор» для переходной характеристики замкнутой
системы исходя из предельно допустимых значений перерегулирования σm и длительности переходного процесса tп (рис. 1.3). При
этом определим длительность переходного процесса как время, по
истечении которого величина курса укладывается в диапазон ±2%
относительно заданного значения Kз.
K/Kз
1+σm
1,02
1
0,98
0
tп
t
Рис. 1.3. Допустимый «коридор» для переходной характеристики
Задача синтеза заключается в нахождении таких значений коэффициентов регулятора kп, kд, kи, чтобы переходная характеристика замкнутой системы укладывалась в допустимый «коридор».
Остальные параметры системы считаются заданными.
8
1.2. Описание программы моделирования
Лабораторная работа выполняется в среде ����������������������
Matlab����������������
6.5 с использованием пакета Simulink.
В соответствии с уравнениями (1.3)–(1.5) формируется Simulinkмодель, показанная на рис. 1.4. Текущим отклонением курса ψ от
заданного значения Kз является выход блока Integrator2, при этом
входной сигнал Kз задается в виде начального условия в блоке Integrator2.
Текущий курс (вход блока Gain2) формируется как разность
между курсом заданным и величиной ψ: K = Kз – ψ. Блок Gain2
служит для нормировки курса к заданному значению с той целью,
чтобы установившееся значение входного сигнала блока NCD�����
��������
����
Outport было равно единице.
В нулевой момент времени текущее значение курса для определенности принимается равным нулю. При
том, что начальные условия для ψ ненулевые (ψ = Kз), система управления будет стремиться свести к нулю значение
ψ = Kз – K, в результате чего величина курса изменяется от нуля до
заданного значения. Тем самым осуществляется маневр по курсу
на величину Kз относительно первоначального значения.
В блоке NCD Outport задается допустимый «коридор» для переходной характеристики. После запуска процедуры настройки параметров регулятора в графической области отображается ход ее
выполнения в виде исходного графика переходного процесса белого
цвета и текущего изменяющегося графика зеленого цвета. В это же
время в командном окне Matlab отображаются текущие значения
настраиваемого параметра и целевой функции оптимизации (функции штрафа). По окончании выполнения процедуры оптимальные
значения настраиваемых переменных, соответствующие кривой
зеленого цвета, сохраняются в рабочем пространстве Matlab.
Отметим, что поскольку в пакете NCD-Blockset используется
градиентный метод поиска минимума функции штрафа, попадание в точку глобального экстремума не гарантируется, даже если
на заданных ограничениях он в принципе достижим. В этом случае
следует изменить начальное приближение для настраиваемого параметра и повторить процедуру настройки.
1.3. Порядок выполнения лабораторной работы
1. Получить вариант задания у преподавателя.
2. Запустить Matlab, создать новый mdl-файл.
9
10
Gain1
b2
Integrator4
1
s
-K-
Saturation
Gain3
Integrator
1
s
Gain7
-K-
MATLAB Fcn1
MATLAB
Function
Integrator1
1
s
Dead Zone
Gain6
kd
Gain4
Ground
Integrator2
1
s
kp
K3
1
s
Gain5
Switch
MATLAB Fcn
MATLAB
Function
Scope
NCD Outport
NCD
OutPort 1
-K-
R2D
Gain2
ki
Radians
to Degrees
Constant
Integrator3
Рис. 1.4. Модель исследуемой системы в пакете Simulink
Gain
b1
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 либо
c3 << c2 );
Initial condition:
-1/sqrt(-c3)
(если судно неустойчиво на курсе и c2 = 0)
− в блок Integrator2
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:
u2 >= Threshold
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
11
Lower saturation limit:
- dmax
в блок Saturation
Upper limit: ddmax
Lower limit: - ddmax
− в блок Gain7
Gain: 1/tau4
− в блок Dead Zone
Start of dead zone:
-dmin
End of dead zone:
dmin
Задаваемые начальные условия соответствуют прямолинейному
движению с установившимся курсом (для устойчивых судов) либо
движению с постоянной скоростью циркуляции (для неустойчивых судов).
4. В меню Simulation выбрать пункт Simulation Parameters и
установить время моделирования 300 с.
5. ������������������������������������������������������������
Создать ����������������������������������������������������
m���������������������������������������������������
-файл в соответствии с прил. 1. Ввести значения параметров в соответствии с вариантом задания (табл. 1.1), обращая
внимание на их размерности.
−
Таблица 1.1
Варианты заданий
Вариант
k1,
с–1
τ1,
с
1
2
3
4
5
6
7
0,03
0,05
0,1
–0,1
–0,12
–0,13
–0,15
10
20
30
–40
–50
–60
–70
τ2, τ3, c2,
с с
с
1
2
3
4
5
6
7
2
5
7
10
10
15
20
c3, с2
δmin,
град
20 20
0,5
35
0
0,3
20 40
0,4
0 –900 0,3
–50 10
0,3
0 –700 0,4
–50 0
0,3
δmax, δ̇ max,
град град/с
35
25
35
35
35
35
25
2
3
4
3
3
4
2
τ4,
с
0,3
0,3
0,25
0,4
0,3
0,3
0,4
Kз, σm, tп,
град % с
60
60
90
90
90
90
90
5
10
10
10
10
20
20
100
90
60
90
90
75
90
6. Вернуться
��������������������������������������������������������
в окно созданной �����������������������������
Simulink���������������������
-модели. Двойным щелчком мыши по блоку NCD Outport раскрыть графическое окно,
представленное на рис. 1.5.
7. Войти в меню Options, выбрать пункт Time Range… и установить значение tmax = 300с.
8. Построить «коридор», в пределах которого должен находится
выходной сигнал исследуемой системы, в соответствии с заданными величинами перерегулирования σm и длительности переходного
процесса tп. Это можно сделать, передвигая красные линии, являющиеся границами «коридора», при помощи мыши. Кроме того,
местоположение этих линий можно установить точно (не в визу12
Рис. 1.5. Графическое окно блока NCD Outport
Рис. 1.6. Построение допустимого «коридора»
для переходной характеристики
альном режиме) при помощи диалоговой панели 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]. Небольшое отрицательное значение координаты
местоположения самой нижней из границ «коридора» задается для
«облегчения» процедуры настройки, так как в противном случае
(при установлении нижней границы точно на нулевом уровне) в мо13
мент времени, равный нулю, переходная характеристика неизбежно касается нижней границы «коридора», в связи с чем функция
штрафа лишается возможности принимать отрицательные значения.
9. В окне созданного mdl-файла войти в меню Simulation�������
������
Param�
eters и установить время моделирования, равное 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.
11. Запустить процедуру настройки коэффициентов регулятора,
нажав виртуальную кнопку Start.
12. Наблюдать развитие процесса настройки коэффициентов по
изменению кривой переходной характеристики (рис. 1.7). Одновременно контролировать ход изменения значения функции штрафа MAX{g} (максимального выхода кривой переходного процесса
за границы «коридора»), отображаемой в командном окне Matlab.
Рис. 1.7. Отображение хода процесса определения
настраиваемых переменных
14
13. Убедиться
�����������������������������������������������������
в завершении процедуры настройки по появлению сообщения “Optimization Converged Successfully” в командном
окне Matlab.
14. Ввести в командном окне Matlab
>> [kp ki kd]
и зафиксировать полученные значения настраиваемых переменных.
15. В случае, если значение MAX{g}, соответствующее последней
итерации, оказалось отрицательным (или равным нулю), считать
задачу синтеза разрешенной, а полученные значения настраиваемых переменных – значениями коэффициентов регулятора kп, kи,
kд. В противном случае откорректировать начальные приближения
(как правило, следует на 10–20% увеличить по модулю величину
kд по отношению к полученному значению), введя новое приближение в командной строке Matlab (kd = …), и повторить процедуру
синтеза, начиная с п. 11.
16. В окне mdl-файла нажать кнопку Start Simulation. Двойным щелчком мыши по блоку Scope раскрыть окно со строящимся
графиком переходного процесса. Изменив масштаб графика, определить амплитуду рыскания в режиме стабилизации.
1.4. Оформление отчета
Отчет о лабораторных исследованиях должен содержать:
− структурную схему исследуемой системы;
− модель системы в пакете Simulink;
− исходные данные для синтеза;
− полученные значения коэффициентов регулятора и соответствующий им график переходного процесса;
− полученное значение амплитуды рыскания в режиме стабилизации.
1.5. Контрольные вопросы
1. Поясните
�������������������������������������������������������
суть процедуры определения настраиваемых переменных с использованием пакета NCD-Blockset.
2. Поясните
�������������������������������������������������������
процедуру линеаризации модели объекта управления.
3. Поясните
������������������������������������������������������
процедуру линеаризации модели рулевого привода.
4. �������������������������������������������������������
Нелинейности какого вида характерны для рулевого привода? Каким образом они проявляются в рассматриваемой системе?
15
5. Запишите
�����������������������������������������������������
передаточную функцию замкнутой линеаризованной системы.
6. Каким образом обеспечивается астатизм в рассматриваемой
системе?
7. Является ли используемый в системе регулятор линейным?
Почему?
8. Чему равно установившееся значение входного сигнала блока
NCD Outport? Почему?
9. ���������������������������������������������������������
Как по графику переходной характеристики определить перерегулирование и длительность переходного процесса?
10. На схеме (рис. 1.4) покажите точки, где можно наблюдать
сигналы x1(t), x2(t), x3(t), u(t), β(t).
16
2. ИССЛЕДОВАНИЕ ТОЧНОСТИ ОПРЕДЕЛЕНИЯ УГЛОВ
ОРИЕНТАЦИИ ПОДВИЖНОГО ОБЪЕКТА
С ИСПОЛЬЗОВАНИЕМ БЕСПЛАТФОРМЕННОЙ
ИНЕРЦИАЛЬНОЙ ВЕРТИКАЛИ
Цель работы: ознакомление с алгоритмом работы бесплатформенной инерциальной вертикали.
2.1. Методические указания по подготовке к работе
Бесплатформенная инерциальная вертикаль предназначена для
выработки углов ориентации маневренных объектов (МО). Основой
вертикали является измерительный блок, содержащий три линейных акселерометра, измеряющих составляющие кажущегося ускорения, и три датчика угловой скорости (ДУС), измеряющих составляющие угловой скорости вращения. Выходные сигналы датчиков
поступают непосредственно в устройство оценивания (УО), которое
определяет мгновенное направление осей чувствительности акселерометров в опорной системе координат. Укрупненная структурная
схема системы представлена на рис. 2.1. Штриховыми линиями на
схеме обозначены связи, которые далее будут искусственно введены для получения математической модели, удобной при синтезе
устройства оценивания.
При построении вычислительного алгоритма бесплатформенной
системы ориентации примем за основу уравнения Пуассона [12,
13], описывающие изменение матрицы направляющих косинусов
cos ψ cos ϑ
sin ϑ
− sin ψ cos ϑ




C =  sin ψ sin γ − cos ψ sin ϑ cos γ cos ϑ cos γ cos ψ sin γ + sin ψ sin ϑ cos γ ,
 sin ψ cos γ + cos ψ sin ϑ sin γ − cos ϑ sin γ cos ψ cos γ − sin ψ sin ϑ sin γ 


где ϑ – угол тангажа; ψ – угол рыскания; γ – угол крена.
v
ξ 4,ξ5,ξ 6
w, ψ, ϑ, γ
МО
ω
Акс
ДУС
wп
ωп
УО
Ƀ γɃ
Ƀ ϑ,
ψ,
ξ1, ξ2,ξ3
Рис. 2.1. Структурная схема инерциальной вертикали
17
В матричной форме эти уравнения имеют вид
E$ 5
 $ 5
EU
(2.1)
где ω – кососимметрическая матрица
    [

  Z
 [
Y
Z 
 Y 


составленная из проекций вектора угловой скорости объекта на
оси, связанной с ним системы координат.
Выражение (2.1) в скалярной форме представляет собой следующую систему дифференциальных уравнений:
C11 = C21ωz − C31ωy , (2.2)
C21 = C31ωx − C11ωz , (2.3)
C31 = C11ωy − C21ωx , (2.4)
C12 = C22 ωz − C32 ωy , (2.5)
C22 = C32 ωx − C12 ωz , (2.6)
C32 = C12 ωy − C22 ωx , (2.7)
C13 = 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) с начальными условиями
18
C12 (0) = sin ϑ0 ,
C22 (0) = cos ϑ0 cos γ 0 ,
C32 (0) = − cos ϑ0 sin γ 0 ,
C13 (0) = − sin ψ 0 cos ϑ0 ,
C23 (0) = cos ψ 0 sin γ 0 + sin ψ 0 sin ϑ0 cos γ 0 ,
C33 (0) = cos ψ 0 cos γ 0 − sin ψ 0 sin ϑ0 sin γ 0 .
Здесь ϑ0, ψ0, γ0 – начальные углы ориентации.
Недостающие элементы матрицы направляющих косинусов
можно определить из соотношений
C11 = C22 C33 − C23 C32 , (2.11)
C21 = C13 C32 − C12 C33 , (2.12)
C31 = C12 C23 − C13 C22 . (2.13)
Тогда параметры ориентации вычисляются по формулам
ϑ = arctg
C12
2
C22
2
+ C32
, ψ = − arctg
C13
C
, γ = − arctg 32 .
C11
C22
Из формул (2.2)–(2.10) следует, что для вычисления параметров
ориентации необходимо знание составляющих угловой скорости
вращения объекта. Такая информация поступает с выходов ДУС с
точностью до ошибок измерений:
п
др
ωx = ωxп + ωxдр , ωy = ωпy + ωдр
y , ωz = ωz + ωz ,
где ωxп, ωyп, ωzп – измеренные (приборные) значения составляющих
угловой скорости вращения, являющиеся входами исследуемой
динамической системы; ωxдр, ωyдр, ωzдр – угловые скорости дрейфа
ДУС. Методические ошибки ДУС будем считать здесь пренебрежимо малыми.
Таким образом, система уравнений относительно направляющих косинусов Cij примет следующий вид:
C12 = C22 ωпz + C22 ωzдр − C32 ωпy − C32 ωдр
(2.14)
y ,
C22 = C32 ωпx + C32 ωxдр − C12 ωzп − C12 ωzдр , (2.15)
19
C32 = C12 ωпy + C12 ωyдр − C22 ωxп − C22 ωxдр , (2.16)
C13 = C23 ωпz + C23 ωzдр − C33 ωпy − C33 ωдр
y ,
(2.17)
C23 = C33 ωпx + C33 ωxдр − C13 ωzп − C13 ωzдр , (2.18)
C33 = C13 ωпy + C13 ωyдр − C23 ωxп − C23 ωxдр . (2.19)
Угловые скорости дрейфа ωxдр, ωyдр, ωzдр будем считать центрированными экспоненциально коррелированными случайными процессами со среднеквадратическими отклонениями σ1, σ2, σ3 и постоянными времени корреляции T1, T2, T3 соответственно. Уравнения
формирующих фильтров для таких процессов имеют вид [3, 14]
 xдр = −T1−1ωxдр + 2T1−1 σ1ξ1, ω
(2.20)
−1 др
−1
 др
ω
y = −T2 ωy + 2T2 σ2 ξ2 , (2.21)
 zдр = −T3−1ωzдр + 2T3−1 σ3 ξ3 , ω
(2.22)
где ξi – некоррелированные друг с другом центрированные белые
шумы единичной интенсивности.
Составляющие абсолютного ускорения объекта wgx, wgy, wgz также представим в виде центрированных экспоненциально коррелированных случайных процессов со среднеквадратическими отклонениями σgx, σgy, σgz и постоянными времени корреляции Tgx, Tgy,
Tgz соответственно:
−1
−1
w gx = −Tgx
wgx + 2Tgx
σ gx ξ4 , (2.23)
−1
−1
w gy = −Tgy
wgy + 2Tgy
σ gy ξ5 , (2.24)
−1
−1
w gz = −Tgz
wgz + 2Tgz
σ gz ξ6 . (2.25)
Выходные сигналы акселерометров wxп, wyп, wzп будут представлять собой зашумленные измерения кажущихся ускорений в проекциях на их оси чувствительности. При движении объекта с постоянной скоростью
20
w xп = w gxC11 + (w gy + g )C12 + w gzC13 + υ1, (2.26)
w yп = w gxC21 + (w gy + g )C22 + w gzC23 + υ 2, (2.27)
п
w z = w gxC31 + (w gy + g )C32 + w gzC33 + υ 3, (2.28)
где uj – некоррелированные друг с другом и с ξi центрированные
белые шумы с заданными интенсивностями Rj; g – ускорение силы
тяжести.
Дифференциальные уравнения (2.14)–(2.25) и алгебраические
уравнения (2.26)–(2.28) представляют собой модель динамической
системы, которая может быть записана в векторно-матричной форме:
(2.29)
Y  G Y V
 # z = h(x) + v, (2.30)
где x = (C12 C22 C32 C13 C23 C33 ωxдр ωyдр ωzдр wgx wgy wgz)T – вектор����
����������
со���
стояния; z = (wxп wyп wzп)T – вектор измерений; u = (ωxп ωyп ωzп)T – ����
вектор входных сигналов; ξ = (ξ1 ξ2 ξ3 ξ4 ξ5 ξ6)T – вектор возмущений;
v = (u1 u2 u3)T – вектор шумов измерений; B ­– матрица возмущений
размерности 12×3, ненулевые элементы которой равны
B7,1 = 2T1−1 σ1, B8,2 = 2T2−1 σ2 , B9,3 = 2T3−1 σ3 ,
−1
−1
−1
B10,4 = 2Tgx
σ gx , B11,5 = 2Tgy
σ gy , B12,6 = 2Tgz
σ gz .
Компоненты вектор-функции f(x,u) будут иметь вид
f1 (x, u) = x2 x9 − x3 x8 + x2u3 − x3u2 ,
f2 (x, u) = x3 x7 − x1x9 + x3u1 − x1u3 ,
f3 (x, u) = x1x8 − x2 x7 + x1u2 − x2u1,
f4 (x, u) = x5 x9 − x6 x8 + x5u3 − x6 u2 ,
f5 (x, u) = x6 x7 − x4 x9 + x6 u1 − x4 u3 ,
f6 (x, u) = x4 x8 − x5 x7 + x4 u2 − x5u1,
f7 (x, u) = −T1−1x7 ,
f8 (x, u) = −T2−1x8 ,
f9 (x, u) = −T3−1x9 ,
21
−1
f10 (x, u) = −Tgx
x10 ,
−1
f11 (x, u) = −Tgy
x11,
−1
f12 (x, u) = −Tgz
x12 .
Компоненты вектор-функции h(x) будут иметь вид
h1 (x) = x10 (x2 x6 − x3 x5 ) + (x11 + g )x1 + x12 x4 ,
h2 (x) = x10 (x3 x4 − x1x6 ) + (x11 + g )x2 + x12 x5 ,
h3 (x) = x10 (x1x5 − x2 x4 ) + (x11 + g )x3 + x12 x6 .
Очевидно, функции f(x,u), h(x) являются нелинейными, так как
их компоненты содержат произведения переменных состояния.
Поэтому для решения задачи фильтрации (оценивания компонент
вектора состояния) будем использовать обобщенный нелинейный
фильтр Калмана-Бьюси. Уравнения фильтра записываются в виде
[14, 15]
x̂ = f (xˆ , u) + K(z − h(xˆ )), (2.31)
K=P
∂f (xˆ , u) −1 R ,
∂x T
(2.32)
∂f (xˆ , u)
∂f T (xˆ , u)
∂h T (xˆ ) −1 ∂h(xˆ )
P =
P
+
P
−
P
R
P + BB T, (2.33)
T
T
∂x
∂x
∂x
∂x
где x̂ – текущая оценка вектора состояния; P – матрица ковариаций вектора ошибок   Y  Ye R – матрица интенсивностей шумов
измерений (диагональная матрица, составленная из элементов Rj);
K – матричный коэффициент усиления фильтра.
∂f (xˆ , u)
Здесь
представляет собой матрицу вида
∂x T
∂f1
∂f1 
 ∂f1
...
 ∂x
∂x12 
 1 ∂x2
 ∂f2
∂f2
∂f2 
...
∂f (x, u) 

x
x
x12 ,
∂
∂
∂
=
2
 1
T
∂x


...


 ∂f12 ∂f12 ... ∂f12 
 ∂x
∂x12 
 1 ∂x2
22
∂f T (xˆ , u)
– транспонированную ма∂x
∂h(xˆ ) ∂h T (xˆ )
трицу. Аналогично определяются производные
,
.
∂x
∂x T
В данном случае
вычисленную в точке x = x̂, а
∂f (x, u)
=
∂x T

x9 +u3−x8 −u2 0
0
0
0
0 −x3 x2 0 0 0 


x3 0 −x1 0 0 0 
0 x7 +u1 0
0
0
 −x9 −u3
 x8 +u2 −x7 −u1 0
0
0
0
0 −x2 x1 0 0 0 0 


0
0
0
0 x9 +u3 −x8 −u2 0 −x6 x5 0 0 0 


0
0
0 −x9 −u3 0 x7 +u1 x6 0 −x4 0 0 0 

0
0
0 x8 +u2 −x7 −u1 0 −x5 x4 0 0 0 0 
=
,
0
0
0
0
0
0 −T1−1 0 0 0 0 0 


0
0
0
0
0
0
0 −T2−1 0 0 0 0 


0
0
0
0
0
0
0 0 −T3−1 0 0 0 

−1

0
0
0
0
0
0
0 0 0 −Tgx
0 0 


−1
0
0
0
0
0
0
0 0 0 0 −Tgy
0 


0
0
0
0
0
0
0 0 0 0 0 −Tgz−1 

∂h(x)
∂x T
=
 x11 + g x6 x10 −x5x10 x12 −x3x10 x2x10 0 0 0 x2x6 x1 x4 


=  −x6 x10 x11 + g x4 x10 x3x10
x12 −x1x10 0 0 0 x3x4 x2 x5 .

x12 0 0 0 x1x5 x3 x6 
 x5x10 −x4 x10 x11 + g −x2x10 x1x10
Начальными условиями для интегрирования уравнений (2.31)–
(2.33) являются априорная оценка вектора состояния x̂(0) и априорная матрица ковариаций вектора ошибок P(0) = P0, где P0 – заданная априорная матрица ковариаций вектора состояния (диагональная матрица, состоящая из дисперсий компонент вектора состояния σ20i ). Оптимальной априорной оценкой вектора состояния
является его безусловное математическое ожидание, которое при
малых углах ϑ0, ψ0, γ0 будем считать приближенно равным вектору xˆ (0) = (0 1 0 0 0 1 0 0 0 0 0 0 0) T. Диагональные элементы
априорной матрицы ковариаций P0 определим как
σ201 = σ202 = σ203 = σ204 = σ205 = 10−4 ;
23
σ206 = 10−6 ; σ207 = σ12 ; σ208 = σ22 ; σ209 = σ23 ;
σ20 10 = σ2gx ; σ20 11 = σ2gy ; σ20 12 = σ2gz .
2.2. Описание программы моделирования
Лабораторная работа выполняется в среде ����������������������
Matlab����������������
6.5 с использованием пакета Simulink.
Для анализа точности оценивания углов ориентации формируется Simulink-модель в соответствии с рис. 2.2. Модели блоков Ob-
MATLAB
Function
angles
x
u
z
em
e(C12),+/-3sigma
Object
em
e(C22),+/-3sigma
e(C32),+/-3sigma
u
x^
z
P
e(C13),+/-3sigma
Kalman Filter
e(C23),+/-3sigma
MATLAB
Function
sigma
MATLAB
Function
angles est
e(C33),+/-3sigma
Scope1
e(omegdx),+/-3sigma
em
theta, deg
e(omegdy),+/-3sigma
e(omegdz),+/-3sigma
em
e(wgx),+/-3sigma
psi, deg
e(wgy),+/-3sigma
e(wgz),+/-3sigma
gamma, deg
Scope
Рис. 2.2. Главное окно модели
24
Scope2
ject (объект) и Kalman Filter (фильтр Калмана), сформированные
как подсистемы (блоки типа Subsystem) в соответствии с уравнениями (2.14)–(2.25) и (2.31)–(2.33), показаны на рис. 2.3 и 2.4 соответственно.
Исходные данные и начальные условия для интегрирования
систем дифференциальных уравнений размещаются в m-файле.
Возмущающие входные воздействия моделируются широкополосными процессами с малой постоянной времени (0,01 с).
Движение объекта задается посредством формирования составляющих угловой скорости, изменяющихся по пилообразному закону с заданными амплитудами Ωx, Ωy, Ωz, частотами Ωx, Ωy, Ωz и
случайной начальной фазой.
K*u
Band-Limited
White Noise
Gain3
1
s
Sine Wave
Sign
-K-
Integrator
1
s
-K-
Integrator1
Sine Wave1 Sign1
1
s
Sine Wave2 Sign2
u(t)
u(t)
Gain1
-K-
Integrator2
MATLAB
Function
Gain
Gain2
omega_x
omega_y
omega_z
2
u
MATLAB
Function
MATLAB Fcn
f(x,u)
dx(t)/dt
x(t)
1
s
Integrator3
MATLAB
Function
h(x)
1
x
3
z
Band-Limited
White Noise1
Рис. 2.3. Блок Object
25
26
Constant
R^-1
2
z
Product
Matrix
Multiply
x^
dh(x^)/dx'
h(x^)
P
P
Constant1
B*B'
Reshape1
Reshape
Integrator1
1
s
dP/dt
Reshape
Reshape
Рис. 2.4. Блок Kalman Filter
MATLAB
Function
1
x^
df(x^,u)/dx'
MATLAB
Function
MATLAB
Function
f(x^,u)
MATLAB
Function
Integrator
1
s
1
u
Math
Function1
uT
Math
Function
uT
2
P
Product1
Matrix
Multiply
Product3
Matrix
Multiply
Product2
Matrix
Multiply
В работе для всех вариантов приняты значения констант
R1 = R2 = R3 = 10–7 м2/с3; T1 = T2 = T3 = 10 с; Tgx = Tgy = Tgz = 1 с;
период изменения составляющих угловой скорости движения объекта 10 с (Ωx = Ωy = Ωz = π/5 с–1).
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))]
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
27
вкладка 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:
В блоках Object и Kalman Filter задать необходимое количество
входов и выходов в соответствии с рис. 2.2.
4. Раскрыть блок Object. Сформировать Simulink���������������
�����������������������
-модель в соответствии с рис. 2.3. Ввести:
− в блок Band-Limited White Noise
Noise power: 1
Sample time: 0.01
− в блок Sine Wave
Amplitude:
1
Bias: 0
Frequency (rad/sec):Omegx
Phase (rad): fi
Sample time: 0
− в блок Sine Wave1
Amplitude:
1
Bias: 0
Frequency (rad/sec):Omegy
Phase (rad): fi
Sample time: 0
− в блок Sine Wave2
Amplitude:
1
Bias: 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
− в блок Integrator2
Initial condition:
(-0.5+abs(fi/pi))*pi/Omegz
− в блок Gain
Gain:Omeg_x*2*Omegx/pi
28
−
в блок 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
− в блок 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
29
−
−
−
−
−
30
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)]
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
Constant value:
B*B’
6. В меню Simulation выбрать пункт Simulation Parameters и
установить время моделирования 50 с.
7. Создать
������������������������������������������������������������
m���������������������������������������������������
����������������������������������������������������
-файл в соответствии с прил. 2. Ввести значения параметров в соответствии с вариантом задания (табл. 2.1), обращая
внимание на их размерности.
Таблица 2.1
Варианты заданий
Вариант
σ1, σ2, σ3, σgx, σgy, σgz, Ωx,
Ωy,
Ωz,
град/с
м/с2
град/с град/с град/с
ϑ0,
град
ψ0,
град
γ0,
град
1
0,01
0,001
6
8
10
1
–1
1
2
0,02
0,002
5
7
12
0
1
–1
3
0,01
0,005
5
10
14
–1
–1
–1
4
0,005
0,001
8
6
16
–2
1
0
5
0,005
0,002
8
5
18
1
–1
–1
6
0,005
0,005
10
15
20
–1
1
–1
7
0,02
0,001
8
10
10
–1
0
1
8. Нажатием клавиши [F5] запустить m-файл на выполнение.
9. ��������������������������������������������������������
Вернуться в главное окно Simulink�����������������������
�������������������������������
-модели. Нажатием кнопки Start Simulation запустить mdl-файл на выполнение.
10. Двойным щелчком мыши по блокам Scope1, Scope2, Scope3
раскрыть окна со строящимися графиками:
31
− углов ориентации ϑ, ψ, γ и их оценок;
− ошибок оценок переменных состояния с оценками точности
(±3σ).
11. Сравнить
������������������������������������������������������
оценки углов ориентации с их истинными значениями. Записать величины максимальных ошибок для отчета.
12. ��������������������������������������������������������
Убедиться в соответствии ошибок оценок переменных состояния оценкам их точности (проконтролировать попадание значений ошибок в коридор ±3σ). Сделать вывод о том, каким образом
ошибки зависят от ориентации объекта.
13. В блоках Band-Limited White Noise, Band-Limited White
Noise1 установить интенсивности входных шумов, равные нулю.
Повторить процедуру моделирования. Определить переменные состояния, переходные процессы для которых имеют затухающий
характер. По длительности переходных процессов оценить время
готовности бесплатформенной вертикали к работе.
2.4. Оформление отчета
Отчет о лабораторных исследованиях должен содержать:
− структурную схему исследуемой системы;
− систему дифференциальных уравнений, описывающих поведение объекта;
− полученные значения экспериментальных данных;
− выводы по работе.
2.5. Контрольные вопросы
1. Какие
���������������������������������������������������������
датчики включает в себя измерительный блок бесплатформенной вертикали? Какие физические величины они измеряют?
2. Какой вид имеют корреляционные функции составляющих
угловой скорости дрейфа ДУС?
3. ������������������������������������������������������
Какие составляющие содержат выходные сигналы акселерометров?
4. Сколько элементов матрицы направляющих косинусов дают
полную информацию об углах ориентации объекта? Какой порядок
имеет исследуемая динамическая система?
5. Какие сигналы являются входными и какие – выходными по
отношению к устройству оценивания?
6. ���������������������������������������������������������
Какое распределение имеет начальная фаза в законе изменения угловой скорости движения объекта?
7. В каком из блоков модели Simulink определяются начальные
условия для интегрирования дисперсионного уравнения?
32
8. ��������������������������������������������������������
В каких блоках модели ����������������������������������
Simulink��������������������������
осуществляется формирование закона изменения составляющих угловой скорости?
9. Каким
���������������������������������������������������
образом моделируются входные возмущающие воздействия типа белого шума?
10. Является ли исследуемая система полностью наблюдаемой?
Почему?
33
3. СИНТЕЗ РЕГУЛЯТОРА ДЛЯ СИСТЕМЫ УПРАВЛЕНИЯ
УПРУГИМ ОБЪЕКТОМ
Цель работы: изучение вопросов построения линейных регуляторов с демпфирующими свойствами.
3.1. Методические указания по подготовке к работе
Стремление создать максимально легкие управляемые конструкции приводит к проявлению их упругих свойств. Учет этих
свойств принципиально важен при проектировании систем автоматического управления подвижными объектами, подверженными
значительным динамическим нагрузкам и силам сопротивления
окружающей среды. Наличие упругости способствует возникновению колебаний в системе управления на различных резонансных
частотах, в ряде случаев приводящих к потере устойчивости и разрушению конструкции.
Созданию эффективных регуляторов, подавляющих колебательные составляющие движения, препятствуют сложность получения достоверной информации об упругих свойствах нелинейного
объекта и значительное изменение резонансных частот во времени.
Применяющийся при решении таких задач метод конечных элементов оказывается непригодным для аэрокосмических аппаратов
сложной формы, состоящих из сотен и тысяч деталей, даже с учетом
широких возможностей современной вычислительной техники.
Одним из возможных путей преодоления возникающих проблем
является переход от исходной системы уравнений в частных производных, описывающих движение упругого объекта, к линеаризованной системе уравнений большей размерности с дальнейшим ее
упрощением [16]. Далее будем считать, что математическая модель
исследуемого объекта для отдельного режима движения получена
с использованием данной методики и записана в виде передаточной
функции
W (s) =
k(s + µ1 )...(s + µ )(s2 + 2δ1ν1s + ν12 )...(s2 + 2δ p ν p s + ν2p )
(s + α1 )...(s + α λ )(s2 + 2d1β1s + β12 )...(s2 + 2dρβρ s + β2ρ )
(3.1)
с постоянными параметрами k, µi, νi, δi, αi, βi, di. При этом коэффициенты αi задают апериодические составляющие движения объекта, а коэффициенты βi, di – колебательные составляющие, или
моды упругих колебаний. Для реальных объектов порядок числителя передаточной функции меньше порядка знаменателя, т.е. l +
2p < λ + 2ρ.
34
Задача управления заключается в подавлении (демпфировании) мод упругих колебаний и получении замкнутой системы с
приемлемыми динамическими характеристиками. Структурная
схема синтезируемой системы управления представлена на рис.
3.1, где g(t) – задающее воздействие; u(t) – управляющее воздействие; w(t) – внешнее возмущение, приведенное к входу системы;
z (t) – выходной сигнал объекта управления; z(t) – наблюдаемый
(измеряемый) сигнал; u(t) – шум измерений; Wр(s) – передаточная
функция регулятора. Возмущения w(t) и u(t) будем считать некоррелированными центрированными белыми шумами с известными
интенсивностями N1 и N2.
Синтез закона управления сложными динамическими системами удобно производить с использованием аппарата пространства
состояний. С целью перехода от передаточной функции объекта к
его описанию в пространстве состояний, в котором каждая из переменных состояния отвечала бы за ту или иную колебательную
либо апериодическую составляющую движения, предварительно
представим передаточную функцию (3.1) в виде суммы простейших дробей. При отсутствии кратных корней у знаменателя такое
представление имеет вид
W (s) =
r
r1
r
+ 2 + ... + n .
s − p1 s − p2
s − pn
Число слагаемых n = λ + 2ρ определяет порядок системы, или
размерность вектора состояния.
В результате разложения передаточной функции W(s) на элементарные дроби образуются две группы слагаемых: с вещественными и с комплексно-сопряженными полюсами [6]. Рассмотрим
каждую из них.
r
Пусть Wi (s) = i
– элементарная передаточная функция с веs − pi
щественным полюсом pi. Такой передаточной функции (от входа u
к выходу z ̃ (i)) соответствует описание в пространстве состояний
w(t)
g(t)
~
z(t)
W(s)
u (t )
Wp (s )
z (t)
υ(t)
Рис. 3.1. Структурная схема замкнутой системы
35
xi = pi xi + u,
z (i) = ri xi .
Здесь z ̃ (i) – i-я аддитивная составляющая выходного сигнала z .
r
r
Теперь предположим, что функции Wi (s) = i и Wi (s) = i +1
s − pi
s − pi +1
имеют комплексно-сопряженные полюса pi и pi+1 (соответственно,
комплексно-сопряженными будут также коэффициенты ri и ri+1).
Чтобы избавиться от мнимой единицы в рассматриваемых передаточных функциях, объединим эти слагаемые в одно:
 (s) = W (s) + W (s) =
W
i
i
i +1
ri
r
−r p − r p + s(ri + ri +1 )
+ i +1 = i2 i +1 i +1 i
.
s − pi s − pi +1
s − ( pi + pi +1 )s + pi pi +1
Передаточной функции Wi (s) соответствует описание в пространстве состояний [1]
xi = xi +1,
xi +1 = − pi pi +1xi + ( 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) с вещественным полюсом, а во втором – одна из функций Wi (s).
~ (1)( s)
W
~ (2)
W ( s)
~
z
...
g+u+w
~ (λ +ρ )( s)
W
Рис. 3.2. Представление структуры объекта в виде
параллельного соединения звеньев
36
Таким образом, математическую модель свободного движения
упругого объекта (при отсутствии задающего воздействия) с учетом
действующего на входе системы возмущения w(t) и шума измерений u(t) можно представить в векторно-матричной форме следующим образом:
x = Ax + B(u + w),
(3.2)
z = Cx + u.
Здесь матрица A имеет блочно-диагональную структуру
 A1

0
A=


 0
0 ...
0 

A2 ...
0 
,

...

0 ... Aλ+ρ 
причем ее блоки A1, …, Aλ+ρ определяются следующим образом:
Ai = pi, если pi – вещественный полюс;
1
 0

Ai = 
, если pi и pi+1 – комплексно-сопряженp
p
p
p
−
+
i
i +1 
 i i +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:
− если 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.
37
Систему уравнений (3.2) можно переписать несколько в ином
виде, удобном для синтеза оптимального регулятора с использованием пакета Matlab:
Y  "Y  #V  
(3.3)
[  $Y  
где ξ – многомерный белый шум с матрицей интенсивностей
Nξ = BBTN1.
Сформируем управление u(t) как выход линейно-квадратичного
гауссовского регулятора (���������������������������������������
LQG������������������������������������
-регулятора) в соответствии с критерием качества
∞

J = M  ∫ (x T (τ)Qx(τ) + u2 (τ))d τ  → min,  0

(3.4)
где Q – задаваемая матрица весовых коэффициентов.
Оптимальный LQG-регулятор представляет собой комбинацию
фильтра Калмана и LQ�����������������������������������������
�������������������������������������������
-регулятора, синтезированного в соответствии с критерием (3.4) [5, 15] и описываемого уравнением
u = −BHxˆ ,
где H – положительно определенное решение алгебраического уравнения Риккати [4]
AT H + HA − HBBT H + Q = 0.
В пакете Matlab�������������������������������������������
�������������������������������������������������
процедура синтеза ������������������������
LQG���������������������
-регулятора автоматизирована и осуществляется с помощью встроенной функции lqg,
аргументами которой являются матрицы, определяющие модель
системы в пространстве состояний, интенсивности входных шумов
и заданные веса в критерии качества.
Матрица Q в квадратичном критерии, как правило, задается
диагональной, что позволяет придать весовым коэффициентам Qjj
достаточно ясный физический смысл, поскольку в этом случае
x T Qx = Q11x12 + Q22 x22 + ... + Qnn xn2 .
А именно, увеличение того или иного весового коэффициента Qjj приведет к дополнительным ограничениям на возможность
увеличения переменной состояния xi по сравнению с другими переменными состояния и сигналом управления. Тем не менее, в том
случае, когда математическая модель объекта управления получена без опоры на физический смысл переменных состояния, выбор
весовых коэффициентов может вызвать существенные трудности.
38
Избежать этих трудностей позволит полученное описание упругого объекта, в котором каждая из переменных состояния формирует
одну из составляющих движения объекта.
Пусть pi и pi+1 – комплексно-сопряженные полюса передаточной
функции W(s). Тогда переменная состояния xi в первом приближении будет описывать гармонические колебания с частотой ωi = Im pi
и некоторой начальной фазой ϕi:
xi (t) ≈ xi max sin(ωi t + ϕi ),
где xi max – амплитуда колебаний. Соответственно, для переменной
xi+1 max запишем
xi +1 (t) ≈ xi max ωi cos(ωi t + ϕi ) = xi +1max cos(ωi t + ϕi ).
Если колебания на частоте ωi  являются преобладающими, выходной сигнал z  будет иметь вид
z(t) ≈ Ci xi (t) + Ci +1xi +1 (t) ≈ (Ci sin(ωi t + ϕi ) + Ci +1ωi cos(ωi t + ϕi ))xi max =
C

=  i sin(ωi t + ϕi ) + Ci +1 cos(ωi t + ϕi )  xi +1max ,
 ωi

а его максимальное значение составит
zmax ≈ xi +1max
Ci2
ω2i
+ Ci2+1 .
Заметим, что передаточная функция от управления u(t) к переменной состояния xi(t) соответствует колебательному звену. Поэтому управление по переменной xi приведет к подавлению не только
упругих колебаний на частоте ωi, но и низкочастотной составляющей спектра входного сигнала, что в свою очередь вызовет увеличение перерегулирования в системе. Руководствуясь данными соображениями, зададим нулевыми весовые коэффициенты Qii при
переменных xi (а также при переменных, формирующих апериодические составляющие движения) и ориентировочно определим коэффициенты Qi+1 i+1 при производных xi+1 = ẋi, воспользовавшись
методом равных вкладов [17]. Согласно данному методу весовые
коэффициенты в критерии (3.4) следует выбирать таким образом,
чтобы максимально допустимые отклонения переменных состояния вносили в функционал качества одинаковый вклад, равный
вкладу максимально допустимого сигнала управления. Положив
приблизительно равным единице коэффициент передачи регулятора на частоте резонанса, получим
39
Qi +1 i +1 =
2
umax
xi2+1max
≈
2
zmax
xi2+1max
≈
Ci2
ωi2
+ Ci2+1 ≈
Ci2
(Im pi )2
+ Ci2+1. (3.5)
Следует учесть, что значения коэффициентов, вычисленных по
формуле (3.5) являются ориентировочными и в ходе выполнения
работы должны быть скорректированы таким образом, чтобы обеспечить выполнение требований к замкнутой системе по быстродействию (длительность переходного процесса не должна превышать
заданного значения tп) и запасу устойчивости (перерегулирование
не должно превышать заданного значения σ). При этом сигнал управления по абсолютной величине не должен превышать заданного
значения umax.
3.2. Описание программы моделирования
Лабораторная работа выполняется в среде ����������������������
Matlab����������������
6.5 с использованием пакета Simulink.
Расчет регулятора и оценка качества переходного процесса выполняются с использованием m-файла. Передаточная функция
объекта определяется посредством задания ее числителя и знаменателя в виде символьных выражений. После запуска файла на выполнение производится построение переходной и амплитудно-частотной характеристик замкнутой системы.
Для исследования работы системы в условиях возмущений формируется Simulink-модель в соответствии с рис. 3.3.
Возмущающие входные воздействия моделируются широкополосными процессами с постоянной времени, меньшей периода ко-
Band-Limited
White Noise
1
Constant
num(s)
den(s)
Transfer Fcn
Scope
x' = Ax+Bu
y = Cx+Du
Scope1
State-Space
Band-Limited
White Noise1
Рис. 3.3. Модель замкнутой системы в пакете Simulink
40
лебаний наиболее высокочастотной колебательной составляющей
движения объекта. Для всех вариантов интенсивности белых шумов принимаются равными 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
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).
41
Таблица 3.1
Варианты заданий
Вариант
σ, %
tп, с
umax
1
60
50
2
2
3
4
5
6
7
50
30
0
60
50
60
10
5
20
25
50
20
W(s)
2
1
2
3
2
4
1
2
3
4
5
6
7
0,3(s + 0,05)(s + 250)(s2 + 180s + 3 ⋅ 104 )(s2 + 90s + 2 ⋅ 104 )(s2 + 150s + 5 ⋅ 104 )
(s + 0,1)(s2 + 10s + 2 ⋅ 105 )(s2 + 5s + 8 ⋅ 104 )(s2 + 0,5s + 104 )(s2 + 0,05s + 4)
0,15(s + 0,5)(s + 1)(s + 100)(s + 250)(s2 + 500s + 105 )(s2 + 30s + 4 ⋅ 104 )
(s + 1,7)(s2 + 10s + 2 ⋅ 105 )(s2 + 5s + 8 ⋅ 104 )(s2 + 0,1s + 600)(s2 + 1,2s + 4)
0,3(s + 1)(s + 250)(s2 + 500s + 105 )(s2 + 0,5s + 100)(s2 + 30s + 4 ⋅ 104 )
(s + 1,7)(s2 + 10s + 2 ⋅ 105 )(s2 + 5s + 8 ⋅ 104 )(s2 + 0,1s + 600)(s2 + 1,2s + 4)
0,1(s + 0,15)(s + 250)(s2 + 28s + 500)(s2 + 85s + 2 ⋅ 104 )(s2 + 300s + 4 ⋅104 )
(s + 0,07)(s2 + 3s + 3 ⋅ 105 )(s2 + 3s + 104 )(s2 + 0,5s + 100)(s2 + 0,05s + 20)
0,2(s + 0,1)(s + 300)(s2 + 70s + 1,5 ⋅ 104 )(s2 + 150s + 4 ⋅ 104 )(s2 + 400s + 5 ⋅ 104 )
(s + 2)(s2 + 7s + 105 )(s2 + 1,5s + 5 ⋅ 104 )(s2 + 0,6s + 5000)(s2 + 0,1s + 6)
0,2(s + 0,05)(s + 200)(s2 + 140s + 2,2 ⋅ 104 )(s2 + 75s + 8000)(s2 + 30s + 4000)
(s + 0,1)(s2 + s + 105 )(s2 + 2,5s + 9000)(s2 + 0,5s + 850)(s2 + 0,05s + 4)
0,15(s + 0,1)(s + 600)(s2 + 50s + 1000)(s2 + 350s + 4 ⋅ 104 )(s2 + 900s + 5 ⋅ 104 )
(s + 2)(s2 + 5s + 105 )(s2 + 0,2s + 5 ⋅ 104 )(s2 + 0,2s + 5000)(s2 + 2s + 6)
6. Нажатием клавиши [F5] запустить файл на выполнение.
7. Проанализировать
�����������������������������������������������������
построенные в раскрывшихся графических окнах переходную (Step Response) и амплитудно-частотную
(���������������������������������������������������������������
Bode�����������������������������������������������������������
����������������������������������������������������������
Diagram���������������������������������������������������
) характеристики. Выделить одну – две наиболее значимые моды колебаний, приблизительно определить их частоты.
42
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).
14. Сравнить
��������������������������������������������������������
полученные в пп. 11–13 значения показателей качества системы управления с заданными в табл. 3.1.
15. В случае? если синтезированная система по какому-либо
показателю не удовлетворяет заданным требованиям, изменить
значения некоторых весовых коэффициентов (кроме последнего,
равного 1) и повторить действия в соответствии с пп. 10–14. При
корректировке весовых коэффициентов следует руководствоваться
следующими соображениями:
− изменению подлежат, прежде всего, коэффициенты при производных составляющих колебательного движения Qi+1 i+1;
− пропорциональное увеличение весовых коэффициентов при
переменных состояния приводит к улучшению качества переходного процесса, но увеличивает управление;
43
− увеличение весового коэффициента при переменной состояния, формирующей одну из мод колебаний, по сравнению с весовым коэффициентом при переменной состояния, формирующей
другую моду, усиливает подавление первой, но ослабляет подавление второй моды;
− наибольшее влияние на показатели качества замкнутой системы, как правило, оказывает коэффициент Qi+1 i+1 при переменной
состояния xi+1, отвечающей за колебания системы с наиболее низкой из резонансных частот;
− с целью уменьшения длительности переходного процесса можно задавать ненулевые весовые коэффициенты при переменных состояния, отвечающих за апериодическое движение, однако перерегулирование в этом случае может сильно возрасти.
16. Определить передаточную функцию регулятора Wр(s), удовлетворяющего всем заданным требованиям. Для этого вызвать
функцию tf(sys2) в командном окне Matlab. Записать полученное
выражение.
3.4. Оформление отчета
Отчет о лабораторных исследованиях должен содержать:
− структурную схему замкнутой системы;
− исходные данные для синтеза;
− значения весовых коэффициентов, использовавшиеся при оптимизации закона управления;
− полученное выражение для передаточной функции регулятора;
− графики переходной и амплитудно-частотной характеристик
исходной и скорректированной систем;
− выводы по работе.
3.5. Контрольные вопросы
1. Как
�������������������������������������������������������
связаны полюса передаточной функции с частотами составляющих колебательного движения?
2. Как изменится вид разложения передаточной функции на
простейшие дроби, если порядок ее числителя будет равен порядку
знаменателя?
3. Все
��������������������������������������������������������
полюса передаточной функции замкнутой системы вещественные. Обязательно ли переходная характеристика будет иметь
апериодический характер? Почему?
4. Для чего в процессе синтеза закона управления передаточная
функция упругого объекта раскладывается на простейшие дроби?
44
5. Передаточная
�������������������������������������������������������
функция объекта с одним входом и одним выходом имеет m нулей и n полюсов. Каковы будут размерности матриц A, B, C при его описании в пространстве состояний?
6. Что
�������������������������������������������������������
можно сказать о полюсах передаточной функции объекта, если его матрица динамики диагональная?
7. Какова особенность использования квадратичного критерия
оптимальности в стохастических системах?
8. С какой целью в работе используется метод равных вкладов?
В чем его суть?
9. ��������������������������������������������������������
Каким будет результат синтеза оптимального �������������
LQG����������
-регулятора, если все весовые коэффициенты при переменных состояния задать равными нулю?
10. Опишите способы формирования линейных стационарных
объектов в среде Matlab.
45
4. СИНТЕЗ АЛГОРИТМА ЭКСТРАПОЛЯЦИИ
УЗКОПОЛОСНОГО СЛУЧАЙНОГО СИГНАЛА
Цель работы: изучение процедуры экстраполяции случайных
процессов с применением алгоритма идентификации.
4.1. Методические указания по подготовке к работе
Управление движением морских судов происходит в условиях
интенсивных волновых возмущений. Эффективным путем уменьшения влияния морского волнения на качество управления является использование как можно более полной информации об угловом
движении судна. В частности, возможность получения экстраполированных (упрежденных) значений углов качки может способствовать предотвращению аварийной ситуации в работе движителя
судна, когда ходовой винт при определенных углах качки и уровне
волнения может оказаться на некоторое время вне водной среды.
Кроме того, решение задачи экстраполяции (прогнозирования)
качки требуется при работе систем управления посадкой корабельной авиации, когда возникает необходимость определения взаимного расположения палубы корабля и летательного аппарата.
Возможные подходы к решению задачи прогнозирования могут быть основаны на построении оптимальных экстраполяторов с
применением теории винеровской или калмановской фильтрации
[8, 18, 19]. Общий недостаток получаемых решений связан с априорным представлением качки в виде стационарного случайного
процесса с известными характеристиками. В то же время известно,
что параметры качки (дисперсия, преобладающая частота, ширина спектра) зависят от целого ряда факторов, таких как балльность
морского волнения, величина курсового угла к волне, скорость
движения судна, направление ветра и т.д. [20]. Отсутствие точных данных о параметрах качки может привести к значительным
ошибкам в прогнозировании. Для получения достоверных оценок
параметров качки с высоким быстродействием может быть использована идентификационная процедура, основанная на нелинейных алгоритмах обработки измерений. Синтез таких алгоритмов
удобно проводить в предположении о возможности представления
процесса качки судна на ограниченных отрезках времени (порядка
нескольких периодов) в виде гармонических колебаний, происходящих с частотой, близкой к преобладающей [21]. При таком описании предполагается, что отклонения амплитуды и частоты качки
от средних значений являются медленноменяющимися случайны46
ми процессами с корреляционными функциями, близкими к экспоненциальным.
Считая несущественными изменения амплитуды колебаний по
сравнению с изменениями частоты, запишем
θ(t) = A 0 sin ((ω0 + ∆ω(t))t + ϕ0 ), (4.1)
где θ(t) – текущее значение угла качки; A0 – амплитуда качки;
ω0 – среднее значение частоты качки; ∆ω(t) – случайные девиации
частоты; ϕ0 – случайная начальная фаза.
Дважды продифференцировав выражение (4.1) по времени с учетом малости производной ∆ω̇ , получим систему дифференциальных
уравнений относительно переменных состояния x1 = θ и x 2 = θ :
x1 = x2 , (4.2)
x2 = −(ω0 + x3 )2 x1, (4.3)
где x3 = ∆ω. Дополним полученную систему уравнением формирующего фильтра для экспоненциально коррелированного процесса
∆ω(t) [14]:
x3 = −αx3 + σω 2αξ, (4.4)
где α – показатель затухания корреляционной функции процесса ∆ω(t); σω – среднеквадратическое отклонение частоты качки от
среднего значения; ξ(t) – входной белый шум формирующего фильтра.
В векторно-матричной форме уравнения (4.2)–(4.4) будут иметь
вид
x = f (x) + Bξ,

где x = (θ θ ∆ω) T – вектор состояния; B =  0 0 σ ω 2α

рица возмущений.
Компоненты вектор-функции f(x) составят
f1 (x) = x2 , T

 – мат
(4.5)
f2 (x) = −(ω0 + x3 )2 x1, (4.6)
f3 (x) = −αx3 + σω 2αξ. (4.7)
Предположим, что измерению подвергаются угол качки и скорость его изменения. При учете случайных (белошумных) погрешностей измерителей уравнения измерений представляются в виде
47
z1 = x1 + υ1, (4.8)
z 2 = x 2 + υ 2, (4.9)
или в векторно-матричной форме
z = Hx + v,
где z = (z1 z2)T – вектор измерений; v = (u1 u2)T – вектор шумов изме1 0 0
рений, H = 
 – матрица наблюдений.
0 1 0
Поскольку функция f(x) является нелинейной, для решения задачи оценивания компонент вектора состояния будем использовать
нелинейный фильтр Калмана-Бьюси в виде [14]
x̂ = f (xˆ ) + K(z − Hxˆ ), (4.10)
K=P
∂f (xˆ ) −1 R ,
∂x T
(4.11)
∂f (xˆ )
∂f T (xˆ )
P =
P
+
P
− PH TR −1HP + BB T, T
∂x
∂x
(4.12)
где x̂ – текущая оценка вектора состояния; P – матрица ковариаций вектора ошибок   Y  Ye R – диагональная матрица интенсивностей шумов измерений (с заданными компонентами R1, R2); K –
матричный коэффициент усиления фильтра.
∂f (xˆ )
В соответствии с выражениями (4.5)–(4.7) производная
∂x T
будет иметь вид
 ∂f1

 ∂x1
∂f (x)  ∂f2
=
∂x T  ∂x1
 ∂f
 3
 ∂x1
∂f1
∂x2
∂f2
∂x2
∂f3
∂x2
∂f1 

∂x3 
0

∂f2  
2
 =  −(ω0 + x3 )
∂x3  
0
∂f3  

∂x3 


0 −2x1 (ω0 + x3 ) . (4.13)

−α
0

1
0
Начальными условиями для интегрирования уравнений (4.10)–
(4.12) являются априорная оценка вектора состояния xˆ (0) = 0  и
априорная матрица ковариаций вектора ошибок P(0) = P0, представляющая собой диагональную матрицу, содержащую средние
квадраты компонент вектора состояния:
48
 A02 2

P0 =  0

 0

0
A02 (ω20
+ σ2ω )
0
0 

2 0 .

σ2ω 
Для получения значения угла качки, упрежденного на время τ,
введем в рассмотрение линеаризованную систему уравнений (4.2),
(4.3) относительно оценок переменных состояния x̂1, x̂ 2 :
xˆ 1 = xˆ 2, (4.14)
xˆ 2 = −(ω 0 + ∆ωˆ ) 2 xˆ 1, (4.15)
где ∆ωˆ = x̂ 3 – текущая оценка девиации частоты, и воспользуемся
формулой Коши при нулевом входном воздействии:
Ye U  
 U  U
Ye U
0
1

где Φ(t + τ, t) = eAτ – фундаментальная матрица; A = 
  –
2
 −(ω 0 + ∆ωˆ ) 0 
матрица динамики системы (4.14), (4.15). Матричная экспонента
может быть найдена с помощью преобразования Лапласа [3]:


F "  - Q&  "
 Для рассматриваемой системы получим
p

pE − A = 
2
 (ω 0 + ∆ωˆ )
p

 p 2 + (ω + ∆ωˆ ) 2
0
( pE − A) −1 = 
 −(ω 0 + ∆ωˆ ) 2
 2
2
 p + (ω 0 + ∆ωˆ )
−1 
;
p 
1

2
p + (ω 0 + ∆ωˆ ) 
;

p

p 2 + (ω 0 + ∆ωˆ ) 2 
2
1

cos(ω 0 + ∆ωˆ )τ
sin(ω 0 + ∆ωˆ )τ 
e Aτ = 
.
ω 0 + ∆ωˆ


ˆ
ˆ
ˆ
(
)sin(
)
cos(
)
−
ω
+
∆ω
ω
+
∆ω
τ
ω
+
∆ω
τ
0
0
0


Отсюда выражение для экстраполированной оценки угла качки
примет вид
49
θˆ (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.
Band-Limited
White Noise
x^
z
Sine Wave
Sine Wave1
u
P
MATLAB
Function
Extrapolation
Kalman Filter
Scope
Band-Limited
White Noise1
theta, rad
Transport
Delay
theta_error, rad
R2D
theta, deg
Radians
to Degrees
R2D
theta_error, deg
Radians
to Degrees1
Scope1
Рис. 4.1. Модель исследуемой системы в пакете Simulink
50
При описании полезного сигнала используется модель в виде
гармонических колебаний с заданной частотой ω1 и начальной фазой ϕ1:
θ(t) = A 0 sin (ω1t + ϕ1 ).
При этом предполагается, что частота ω1 может существенно отличаться от частоты настройки фильтра ω0.
Возмущающие входные воздействия моделируются широкополосными процессами с малой постоянной времени (0,01 с).
В блоке Scope1 можно пронаблюдать полезный сигнал, его оценку, приведенную ко времени экстраполяции, а также ошибку оценивания. В блоке Scope формируются графики изменения во времени элементов ковариационной матрицы P.
4.3. Порядок выполнения лабораторной работы
1. Получить вариант задания у преподавателя.
2. Запустить Matlab, создать новый mdl-файл.
3. Сформировать Simulink-модель в соответствии с рис. 4.1.
Ввести:
− в блок Sine Wave
Amplitude:
A0
Bias: 0
Frequency (rad/sec): omeg1
Phase (rad): fi1
Sample time: 0
− в блок Sine Wave1
Amplitude:
A0*omeg1
Bias: 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)
51
−
в блок Transport Delay
Time delay: tau
− в блок Scope1 (меню ′Scope′ parameters)
вкладка General:
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
Output dimensions:[3,3]
− в блок Reshape1
Output dimensionality:
Customize
Output dimensions:[2,3]
− в блоки Math Function, Math Function1
Function:
transpose
− в блок Constant1
Constant value:
B*B’
52
5. В меню Simulation выбрать пункт Simulation Parameters и
установить время моделирования 20 с.
6. ������������������������������������������������������������
Создать ����������������������������������������������������
m���������������������������������������������������
-файл в соответствии с прил. 4. Ввести значения параметров в соответствии с вариантом задания (табл. 4.1), обращая
внимание на их размерности. Задать ω1 = 0,5 с–1, τ = 1 с.
Таблица 4.1
Варианты заданий
Вариант A0, град
1
2
3
4
5
6
7
3
4
5
5
6
7
8
ϕ1, град R1, град2⋅c R2, град2/c
0
45
60
90
225
180
270
10–6
10–4
10–5
10–6
10–5
10–6
10–5
10–6
2⋅10–4
5⋅10–5
2⋅10–4
10–4
2⋅10–4
5⋅10–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 раскрыть окно со
строящимися графиками:
− полезного сигнала θ и его экстраполированной оценки θ̂;
− ошибки экстраполяции ε θ = θ − θ̂.
10. Определить максимальную ошибку экстраполяции εm в установившемся режиме.
11. Повторить процедуру моделирования, задавая значения ω1
и τ в соответствии с табл. 4.2. Зафиксировать полученные значения
ошибок экстраполяции. Проанализировать зависимость точности
экстраполяции от частоты полезного сигнала.
Таблица 4.2
Максимальная ошибка экстраполяции εm, град
ω1,
τ, с
с–1
0,5
1
2
3
4
5
6
1
3
5
4.4. Оформление отчета
Отчет о лабораторных исследованиях должен содержать:
− математическую модель динамики системы;
53
− модель реального сигнала;
− заполненную табл. 4.2;
− выводы по работе.
4.5. Контрольные вопросы
1. ��������������������������������������������������������
Каковы достоинства и недостатки алгоритмов прогнозирования, основанных на теории оптимальной фильтрации?
2. Каким
�������������������������������������������������������
образом может быть описана качка судна на нерегулярном волнении?
3. Какие параметры характеризуют качку судна?
4. Для
��������������������������������������������������������
чего в рассмотренном алгоритме экстраполяции используется идентификационная процедура?
5. Почему действительная матрица ковариаций вектора ошибок
может отличаться от расчетной, вырабатываемой в фильтре Калмана?
 + ω20 x = 0? Какие
6. Какой вид имеет общее решение уравнения x
начальные условия необходимо задать, чтобы найти его частное решение?
7. Что такое фундаментальная матрица линейной системы? При
решении каких задач она используется?
8. Как
�����������������������������������������������������
найти фундаментальную матрицу линейной стационарной системы?
9. Укажите причины возрастания ошибок экстраполяции при
увеличении времени упреждения.
10. Какое соотношение амплитуд и какую разность фаз имеют
гармонический сигнал и его производная?
54
Библиографический список
1. Андриевский Б. Р., Фрадков А. Л. Избранные главы теории автоматического управления с примерами на языке MATLAB. СПб.:
Наука, 1999. 467 с.
2. Андриевский Б. Р., Фрадков А. Л. Элементы математического
моделирования в программных средах MATLAB 5 и Scilab. СПб.:
Наука, 2001. 50 с.
3. �������������������������������������������������������
Методы классической и современной теории автоматического управления: Учебник. Т.1: Анализ и статистическая динамика
систем автоматического управления / Под ред. Н. Д. Егупова. М.:
Изд-во МГТУ им. Н. Э. Баумана, 2000. 748 с.
4. Методы классической и современной теории автоматического
управления: Учебник. Т.2: Синтез регуляторов и теория оптимизации систем автоматического управления / Под ред. Н. Д. Егупова.
М.: Изд-во МГТУ им. Н. Э. Баумана, 2000. 736 с.
5. Методы классической и современной теории автоматического
управления: Учебник. Т.3: Методы современной теории автоматического управления / Под ред. Н. Д. Егупова. М.: Изд-во МГТУ им.
Н. Э. Баумана, 2000. 748 с.
6. Сергиенко А. Б. Цифровая обработка сигналов. СПб.: Питер,
2003. 604 с.
7. Ерофеев А. А. Теория автоматического управления: Учебник
для вузов. СПб.: Политехника, 2002. 302 с.
8. Радиоавтоматика:
�����������������������������������������������������������
Учеб. пособие / Под ред. В. А. Бесекерского. М: Высш. школа, 1985. 271 с.
9. Вагущенко Л. Л., Цымбал Н. Н. Системы автоматического управления движением судна. Одесса: ЛАТСТАР, 2002. 310 с.
10. Березин С. Я., Тетюев Б. А. Системы автоматического управления движением судна по курсу. Л.: Судостроение, 1990. 256 с.
11. Соболев Г. В. Управляемость корабля и автоматизация судовождения. Л.: Судостроение, 1976. 480 с.
12. Анучин О. Н., Емельянцев Г. И. Интегрированные системы
ориентации и навигации для морских подвижных объектов //
СПб.: ГНЦ РФ – ЦНИИ “Электроприбор”, 2003. 390 с.
13. Бромберг П. В. Теория инерциальных систем навигации. М.:
Наука, 1979. 278 с.
14. Степанов О. А. Применение теории нелинейной фильтрации в задачах обработки навигационной информации // СПб: ГНЦ
РФ – ЦНИИ “Электроприбор”, 1998. 370 с.
15. Брайсон А., Хо Ю-Ши. Прикладная теория оптимального управления: Пер. с англ. / под ред. А. М. Летова. М.: Мир, 1972. 544 с.
55
16. Бродский С. А., Небылов А. В., Панферов А. И. Контроль параметров и стабилизация упругих колебаний сложных объектов и
конструкций // Гироскопия и навигация. 2006. №4(55). С. 105.
17. Колганов А. Р. Основные разделы современной теории автоматического управления: Электронный конспект лекций. – http://
elib.ispu.ru/library/lessons/kolganov2/index.html
18. Ривкин С. С. Прогнозирование положения палубы корабля
на качке // Вопросы кораблестроения. Сер. Навигация и гироскопия. 1980. Вып. 49. С. 55–64.
19. Крайнов В. И., Тупысев В. А. Об экстраполяции вырабатываемых значений углов качки корабля // Гироскопия и навигация. 1994. №1(4). С. 58–64.
20. Бородай И. К., Нецветаев Ю. А. Качка судов на морском волнении. Л.: Судостроение, 1969. 432 с.
21. Максимова Г. Н. Экстраполирование узкополосного сигнала // Изв. вузов. Приборостроение. 1973.
56
ПРИЛОЖЕНИЕ 1
Текст m-файла к лабораторной работе №1
%Начальные приближения для коэффициентов ПИДрегулятора
kpid
kp =
ki =
kd =
= [-1 -0.01 -1];
kpid(1);
kpid(2);
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; %отклонение курса заданного, рад
57
ПРИЛОЖЕНИЕ 2
Текст m-файла к лабораторной работе №2
clear; %очистка рабочей области
close all; %закрытие окон
rand(‘state’,0); %установка датчика случайных чисел в
g = 9.81; %СКО
sig1
sig2
sig3
%исходное состояние
%ускорение силы тяжести
угловой скорости дрейфа, 1/с
= 0.005/180*pi;
= 0.005/180*pi;
= 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);
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/с
58
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);
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);
59
ПРИЛОЖЕНИЕ 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+10
00)*(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+5
E+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 %если полюс комплексный, а предыдущий
%был вещественным
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);
60
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; %новое окно
bodemag(sys1,’r--’); %ЛАХ разомкнутой системы
hold on; %в том же окне
bodemag(sys); %ЛАХ замкнутой системы
61
ПРИЛОЖЕНИЕ 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;
R2 = 5E-5*(pi/180)^2;
%интенсивность шума измерений
%угла качки, рад^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]); %априорная матрица ковариаций
62
Содержание
Предисловие.................................................................... 1. Синтез системы стабилизации судна по курсу с заданным
качеством переходного процесса ......................................... 1.1. Методические указания по подготовке к работе........... 1.2. Описание программы моделирования........................ 1.3. Порядок выполнения лабораторной работы................ 1.4. Оформление отчета................................................. 1.5. Контрольные вопросы............................................. 2. Исследование точности определения углов ориентации
подвижного объекта с использованием бесплатформенной
инерциальной вертикали................................................... 2.1. Методические указания по подготовке к работе........... 2.2. Описание программы моделирования........................ 2.3. Порядок выполнения лабораторной работы................ 2.4. Оформление отчета................................................. 2.5. Контрольные вопросы............................................. 3. Синтез регулятора для системы управления упругим
объектом......................................................................... 3.1. Методические указания по подготовке к работе........... 3.2. Описание программы моделирования........................ 3.3. Порядок выполнения лабораторной работы................ 3.4. Оформление отчета................................................. 3.5. Контрольные вопросы............................................. 4. Синтез алгоритма экстраполяции узкополосного
случайного сигнала........................................................... 4.1. Методические указания по подготовке к работе........... 4.2. Описание программы моделирования........................ 4.3. Порядок выполнения лабораторной работы................ 4.4. Оформление отчета................................................. 4.5. Контрольные вопросы............................................. Библиографический список................................................ Приложение 1.................................................................. Текст m-файла к лабораторной работе №1............................. Приложение 2.................................................................. Текст m-файла к лабораторной работе №2............................. Приложение 3.................................................................. Текст m-файла к лабораторной работе №3............................. Приложение 4.................................................................. Текст m-файла к лабораторной работе №4............................. 3
4
4
9
9
15
15
17
17
24
27
32
32
34
34
40
41
44
44
46
46
50
51
53
54
55
57
57
58
58
60
60
62
62
63
Учебное издание
Панферов Александр Иванович
Лопарев Алексей Валерьевич
КОМПЬЮТЕРНЫЙ АНАЛИЗ И
СИНТЕЗ СИСТЕМ ОРИЕНТАЦИИ,
СТАБИЛИЗАЦИИ И НАВИГАЦИИ
Учебно-методическое пособие
Редактор В. П. Зуева
Верстальщик С. Б. Мацапура
Сдано в набор 10.11.08. Подписано к печати 17.11.08.
Формат 60×84 1/16. Бумага офсетная. Печать офсетная.
Усл.-печ. л. 3,72. Уч.-изд. л. 3,97. Тираж 200 экз. Заказ № 670
Редакционно-издательский центр ГУАП
190000, Санкт-Петербург, Б. Морская ул., 67
Документ
Категория
Без категории
Просмотров
0
Размер файла
1 546 Кб
Теги
panferov, 05b50e71db
1/--страниц
Пожаловаться на содержимое документа