Математическая модель в пространстве состояний линейного стационарного объекта управления

СОДЕРЖАНИЕ

1. Анализ объекта управления

1.1 Анализ линейного стационарного объекта управления, заданного передаточной функцией

1.2 Получение математической модели в пространстве состояний линейного стационарного объекта управления, заданного передаточной функцией

1.2.1 Матрица Фробениуса

1.2.2 Метод параллельной декомпозиции

2. Решение задачи быстродействия симплекс-методом

3. Оптимальная l – проблема моментов

3.1 Оптимальная l – проблема моментов в пространстве “вход-выход”

3.2 Оптимальная l – проблема моментов в пространстве состояний

4. Нахождение оптимального управления с использованием грамиана управляемости (критерий – минимизация энергии)

5. Аналитическое конструирование оптимальных регуляторов (акор)

5.1 Стабилизации объекта управления на полубесконечном интервале времени

5.1.1 Решение алгебраического уравнения Риккати методом диагонализации

5.1.2 Решение алгебраического уравнения Риккати интегрированием в обратном времени до установившегося состояния

5.2 Стабилизации объекта управления на конечном интервале времени

5.3 Задача акор – стабилизации для компенсации известного возмущающего воздействия.

5.4 Задача акор для отслеживания известного задающего воздействия. i подход

5.5 Задача акор для отслеживания известного задающего воздействия. ii подход (линейный сервомеханизм)

5.6 Задача акор – слежения со скользящими интервалами.

6. Синтез наблюдателя полного порядка

Литература

Приложение

PlotTimeFrHaract. m

ProstranstvoSostoyanii. m

SimplexMetod2.m

Optimal_L_problem_moments. m

Gramian_Uprav. m

AKOR_stabilizaciya_na_polybeskon_interval. m

AKOR_stabilizaciya_na_konech_interval. m

Sravnenie_stabilizacii. m

AKOR_stabilizaciya_pri_vozmusheniyah. m

AKOR_slegenie_na_konech_interval_I_podxod. m

AKOR_slegenie_na_konech_interval_II_podxod. m

AKOR_slegenie_so_skolz_intervalami_Modern. m

Sintez_nablyud_polnogo_poryadka. m

Solve_Riccati_Method_Diag. m

Solve_Riccati_Method_Revers_Integr. m

Vozmyshyayushee_Vozdeistvie_Discrete_Revers. m

Zadayushee_Vozdeistvie_Discrete_Revers_Modern. m

1. Анализ объекта управления

1.1Анализ линейного стационарного объекта управления, заданного передаточной функцией

Передаточная функция данного объекта имеет вид:

,

Где:

, ;

, , , , , .

Или

.

Нули передаточной функции:

Полюса передаточной функции (полученные стандартными функциями среды Matlab 7.4):

Рис.1. График расположения нулей и полюсов передаточной функции объекта на комплексной плоскости.

Найдем временные характеристики объекта управления.

К временным характеристикам относятся и .

– переходная характеристика;

– импульсная переходная функция;

Для нахождения и воспользуемся пакетом Matlab 7.4.

,

Аналитическое выражение для :

В этом случае имеет вид

Рис.2. График переходной характеристики .

Рис.3. График переходной характеристики на интервале (увеличенное).

,

Аналитическое выражение для :

.

В этом случае имеет вид

Рис.4. График импульсной переходной характеристики .

Рис.5. График импульсной переходной характеристики На интервале (увеличенное).

Найдем частотные характеристики объекта управления.

К частотным характеристикам относятся:

Амплитудно – частотная характеристика (АЧХ),

Фазо – частотная характеристика (ФЧХ),

Амплитудно – фазовая частотная характеристика (АФЧХ),

Аналитическое выражение для АЧХ:

.

В этом случае АЧХ имеет вид

Рис.6. График АЧХ

Рис.7. График АЧХ на интервале (увеличенное). Аналитическое выражение для ФЧХ:

В этом случае ФЧХ имеет вид

Рис.8. График ФЧХ.

Рис.9. График ФЧХ на интервале (увеличенное).

Рис.10. График АФЧХ.

Рис.11. График АФЧХ (увеличенное).

Аналитическое выражение для ЛАЧХ:

.

В этом случае ЛАЧХ имеет вид

Рис.12. График ЛАЧХ.

Аналитическое выражение для ЛФЧХ:

В этом случае ЛФЧХ имеет вид

Рис.13. График ЛФЧХ.

1.2 Получение математической модели в пространстве состояний линейного стационарного объекта управления, заданного передаточной функцией

Передаточная функция данного объекта имеет вид:

,

Где:

, ;

, , , , , .

Или

Описание системы в пространстве состояний имеет следующий вид:

Переходя в область изображений описание системы в пространстве состояний будет иметь следующий вид:

1.2.1 Матрица Фробениуса

Получим выражения, которые определяют вектор состояний и выход заданного объекта в общем виде:

.

.

Тогда получим:

(1)

(2)

Числитель передаточной функции имеет вид: .

Знаменатель передаточной функции:

.

Тогда согласно равенству (1) и (2) имеем

,

.

Перейдем из области изображений в область оригиналов

,

И затем перейдем к нормальной форме Коши

.

Запишем матрицы состояний

, ,

Численное значение матриц состояний:

, ,

1.2.2 Метод параллельной декомпозиции

Запишем передаточную функцию объекта в другом виде, а именно:

Или

.

Согласно формуле получим

Рассмотрим каждое из слагаемых в отдельности согласно принципу параллельной декомпозиции.

A. ,

.

B. ,

.

C. ,

,

,

D. ,

Получим выход системы:

Запишем матрицы состояний

, ,

Вычисление коэффициентов разложения дробной рациональной функции на сумму элементарных дробей и проверка правильности получения матриц состояния сделано с помощью пакета Matlab 7.4 (скрипт ProstranstvoSostoyanii. m)

Получены следующие результаты:Матрица СЛАУ:

, ,

,

Численное значение матриц состояний:

, ,

.

2. Решение задачи быстродействия симплекс-методом

Дана система:

(3)

1. Проверим управляемость данной системы.

Запишем систему ДУ в матричном виде:

,

Где .

Данная система является стационарной, ее порядок , поэтому матрица управляемости имеет вид:

Найдем матрицу управляемости:

Ранг матрицы управляемости равен порядку системы, следовательно, данная система является управляемой.

следовательно .

Собственные числа матрицы найдем из уравнения :

Действительные части собственных значений матрицы являются неположительными, следовательно, все условия управляемости выполнены.

2. Ссылаясь на решение задачи быстродействия из ДЗ№2 по СУЛА “Решение задачи быстродействия” имеем:

Запишем зависимости , , полученные при решении систем дифференциальных уравнений:

:

:

:

:

Перейдем к дискретной модели заданной системы. Имеем

(4)

Где шаг дискретизации и соответствующие матрицы

(5)

Пусть управление ограничено интервальным ограничением

(6)

Тогда на шаге имеем

(7)

Известны начальная и конечная точки

Где – оптимальное число шагов в задаче быстродействия.

Решается задача быстродействия

А) Формирование задачи быстродействия как задачи линейного программирования

Конечная точка в дискретной модели представлена в виде

(8)

Получаем – равенств

(9)

Для приведения ограничений (9) к канонической форме сделаем необходимое преобразование в правой и левой частях, чтобы правые части были неотрицательными (если правая часть меньше нуля, то домножаем на (-1) левую и правую части). Отметим проведенные изменения точкой в правом верхнем углу соответствующих векторов

. (10)

Для того чтобы получить необходимый допустимый базис для задачи линейного программирования, добавим формально остаточные искусственные переменные (). Таким образом, уравнения (10) представляются в виде

(11)

Так как текущее управление – управление имеет любой знак, то сделаем необходимую замену

Тогда уравнения (11) примут вид

(12)

Введем остаточные переменные в ограничения на управление

(13)

При объединении выражений (12) и (13) получаем ограничений.

Начальный допустимый базис состоит из остаточных и остаточных искусственных переменных

Формируем целевую функцию (по второму методу выбора начального допустимого базиса)

(14)

Б) Решение задачи быстродействия

Предположим, что , где – оптимальное число шагов. Так как значение нам неизвестно (но известно точно), выбираем некоторое начальное и решаем задачу линейного программирования (12)-(14).

При этом

Общее число столбцов в симплекс-таблице:

Число базисных переменных:

Сформируем Строку. Имеем

Выразим из уравнения (12) начальные базисные переменные

И подставим в целевую функцию. Получим – строку

(15)

Решаем задачу (12) – (14) симплекс-методом.

В случае,

Если , – малое число

Иначе

1) если Увеличить и целое, рвернуться к первому шагу формирования задачи линейного программирования;

2) если (не все управления будут равны предельным, могут быть, в том числе нулевые)), , уменьшить , вернуться к первому шагу формирования задачи линейного программирования.

Решения данной задачи получено с помощью пакета Matlab 7.4 (скрипт SimplexMetod2.m):

Рис. 14 . График фазовой координаты .

Рис. 15 . График фазовой координаты .

Рис. 16 . График .

Рис. 17 . График оптимального управления .

Выводы: Сравнивая полученные результаты с результатами полученными в ДЗ№2 по СУЛА, можно сделать вывод, что решения совпадают, с точностью до .

3. Оптимальная L – проблема моментов 3.1 Оптимальная L – проблема моментов в пространстве “вход-выход”

Укороченная система данного объекта имеет вид:

,

Где:

;

;

;

;

;

.

Полюса укороченной передаточной функции:

;

;

;

;

.

Заданы начальные и конечные условия:

, , .

Для определения начальных и конечных условий для воспользуемся следующей формулой:

,

Где матрица имеет следующий вид

,

Где , .

ИПФ укороченной системы:

Составим фундаментальную систему решений:

ФСР: .

Составим матрицу .

, где – матрица Вронского

,

Тогда

.

Составим моментные уравнения (связь между входом и выходом):

Моментные функции определяются по следующей формуле

Составим моментные функции:

Найдем моменты по следующей формуле:

.

Числовое значение найденных моментов:

Составим функционал качества, который имеет следующий вид:

При условии, что :, т. е.

Выразим из данного условия , тогда получим следующее равенство:

.

Подставляя полученное равенство в функционал и заменяя их правыми частями получаем

Найдем частные производные и приравняем их к нулю. Решая полученную систему уравнений, определяем оптимальные значения коэффициентов , а вычислим по формуле

.

Т. о. имеем:

Минимальная энергия:

Найдем управление по следующей формуле:

Тогда оптимальное управление

.

3.2 Оптимальная L – проблема моментов в пространстве состояний

Система задана в виде:

Решение ДУ имеет вид:

, при имеем:

.

Составим моментные уравнения:

Подставляя необходимые данные в выше приведенные формулы, получим следующие моменты и моментные функции:

Числовое значение найденных моментов:

Моментные функции:

Заметим, что моменты и моментные функции совпадают с моментами и моментными функциями, найденными в пункте (а).

Из этого следует, что функционал, значения , управление и минимальная энергия будут иметь точно такие же числовые значения и аналитические выражения, как и в пункте (3.1).

Оптимальное управление имеет вид:

Проверим правильность полученного решения.

Эталонные значения координат в начальный и конечный момент времени:

,

,

Найденные значения координат в начальный и конечный момент времени:

,

,

Вычислим погрешность полученных результатов:

,

,

Ниже представлены графики полученного решения с помощью скрипта Optimal_L_problem_moments. m.

Рис. 18. Графики фазовых координат системы при переходе из в .

Рис. 19. Графики выходных координат системы при переходе из в .

Рис.20. График оптимального управления .

Выводы: Задача перевода системы из начальной точки в конечную с помощью L-проблемы моментов в пространстве состояний и в пространстве вход-выход была решена с точностью до 12-го знака после запятой. Результаты, полученные при переводе системы из начальной точки в конечную, полностью совпадают.

4. Нахождение оптимального управления с использованием грамиана управляемости (критерий – минимизация энергии)

Система имеет вид:

С начальными условиями:

,

.

Составим матрицу управляемости и проверим управляемость системы:

.

Составим грамиан управляемости для данной системы:

Найдем грамиан по формуле:

Тогда управление имеет вид:

.

Или

Ниже представлен график оптимального управления полученного с помощью скрипта Gramian_Uprav. m.:

Рис.21. График оптимального управления .

Графики фазовых координат аналогичны, как и в оптимальной L – проблеме моментов.

Сравним управление, полученное в начальной и конечной точках в пунктах 3 и 4 соответственно:

и

Выводы: Как видно, значения граничных управлений совпадают. А это значит, что задача перевода объекта из начального состояния в конечное решена с высокой степенью точности и с минимальной энергией.

Графическое сравнение оптимальных управлений из пунктов 3 и 4:

Рис.21. Сравнение графиков оптимального управления .

5. Аналитическое конструирование оптимальных регуляторов (АКОР) 5.1 Стабилизации объекта управления на полубесконечном интервале времени

Рассмотрим линейный объект управления, описываемый системой дифференциальных уравнений в нормальной форме

Необходимо получить закон управления

Минимизирующий функционал вида

Начальные условия для заданной системы

Моменты времени фиксированы. Матрицы – симметричные неотрицательно определенные:

Матрица – положительно определенная:

Матричное дифференциальное уравнение Риккати имеет вид:

Если линейная стационарная система является полностью управляемой и наблюдаемой, то решение уравнения Риккати при стремится к установившемуся решению не зависящему от и определяется следующим алгебраическим уравнением:

В рассматриваемом случае весовые матрицы и в функционале не зависят от времени.

Оптимальное значение функционала равно

И является квадратичной функцией от начальных значений отклонения вектора состояния.

Таким образом, получаем, что при оптимальное управление приобретает форму стационарной обратной связи по состоянию

Где – решение алгебраического матричного уравнения Риккати.

5.1.1. Решение алгебраического уравнения Риккати методом диагонализации

Для решения данной задачи найдем весовые матрицы и :

Выберем произвольно , тогда

Взяв значения из решения задачи L – проблемы моментов получим:

Матрицы системы имеют вид:

, .

Введем расширенный вектор состояния .

Тогда матрица Z будет иметь следующий вид: ,

Или в численном виде

.

Собственные значения матрицы : .

Зная собственные значения и собственные вектора матрицы Z, построим матрицу

По определению все решения должны быть устойчивы при любых начальных условиях , т. е. при . Чтобы не оперировать комплексными числами, осуществим следующий переход. Пусть:

Тогда матрица формируется следующим образом:

.

Можно показать, что матрицу можно получить из прямой матрицы собственных векторов:

,

.

Установившееся решение уравнения Риккати, полученное с помощью скрипта Solve_Riccati_Method_Diag. m. имеет вид:

5.1.2 Решение алгебраического уравнения Риккати интегрированием в обратном времени до установившегося состояния

Весовые матрицы И Такие же как и в пункте (5.1.1).

Матрицы тоже аналогичны.

Запишем уравнение Риккати

.

Зная, что , решаем уравнение методом обратного интегрирования на достаточно большом интервале (примерно 10 с.), получим установившееся решение с помощью скрипта

Solve_Riccati_Method_Revers_Integr. m.:

Рис.22. Графики решения уравнения Риккати.

Найдем разницу между решениями уравнения Риккати в пунктах 5.1.1 и 5.1.2:

Выводы: сравнивая решения полученные в пунктах 5.1.1 и 5.1.2 можно сказать, что решения уравнения Риккати первым и вторым методами совпадают с заданной точностью. Погрешность расхождения решений невелика.

Используя скрипт AKOR_stabilizaciya_na_polybeskon_interval. m получим коэффициенты регулятора, фазовые координаты системы и управление.

Рис.23. Графики коэффициентов регулятора обратной связи.

Рис.24. Графики фазовых координат.

Рис.25. График управления.

Выводы: т. к. решения уравнения Риккати методом диагонализации и интегрирования в обратном времени дают практически одинаковый результат, то можно считать, что задача АКОР – стабилизации на полубесконечном интервале решена с заданной точностью.

5.2 Стабилизации объекта управления на конечном интервале времени

Рассмотрим линейный объект управления, описываемый системой дифференциальных уравнений в нормальной форме

Начальные условия для заданной системы

Время стабилизации .

Необходимо получить закон управления

Минимизирующий функционал вида

Закон оптимального управления в данной задаче имеет вид

Матричное дифференциальное уравнение Риккати будет иметь следующий вид:

Если обозначить то можно записать

Уравнение замкнутой скорректированной системы примет вид

Матрицы Заданы в пункте 5.1.1.

Весовые матрицы И Имеют следующий вид:

, .

Используя скрипт AKOR_stabilizaciya_na_konech_interval. m получили следующие результаты:

Рис.26. Графики решения уравнения Риккати.

Рис.27. Графики коэффициентов регулятора обратной связи.

Рис.28. Графики фазовых координат.

Рис.29. График управления.

Сравним, как стабилизируется система управления с постоянными и переменными коэффициентами регулятора обратной связи на начальном этапе:

Рис.30. Графики фазовых координат.

Выводы: из графиков видно, что система, у которой коэффициенты регулятора меняются со временем, стабилизируется не хуже, чем, система, у которой коэффициенты регулятора не изменяются.

5.3 Задача АКОР – стабилизации для компенсации известного возмущающего воздействия

Рассмотрим систему вида

,

Где – возмущающее воздействие.

Матрицы Заданы в пункте 5.1.1.

Весовые матрицы И Имеют следующий вид:

, .

Начальные условия для заданной системы .

Время стабилизации .

Задаем возмущающее воздействие только на первую координату, так как только она имеет значение

и .

Решение задачи стабилизации сводится к решению уравнения Риккати

С начальными условиями:

Введем вспомогательную вектор-функцию , ДУ которой имеет вид:

С начальными условиями: .

Управление определяется по формуле:

.

Используя скрипт AKOR_stabilizaciya_pri_vozmusheniyah. m, получили следующие результаты:

Рис.31. Графики решения уравнения Риккати.

Рис.32. Графики коэффициентов регулятора обратной и прямой связи.

Рис.33. График возмущающего воздействия.

Рис.34. График вспомогательной вектор – функции.

Рис.35. Графики фазовых координат.

Рис.36. График управления.

Рис.37. График возмущающего воздействия.

Рис.38. График вспомогательной вектор – функции.

Рис.39. Графики фазовых координат.

Рис.40. График управления.

Выводы: По графикам фазовых координат при различных воздействиях видно, что влияние возмущающего воздействия не существенно и фазовые координаты устанавливаются в ноль. При этом видно, что графики первой фазовой координаты при различных воздействиях мало отличаются друг от друга, т. е. система отрабатывает любое возмущение.

5.4 Задача АКОР для отслеживания известного задающего воздействия. I подход

Система задана в виде:

Матрицы Заданы в пункте 5.1.1.

Весовые матрицы И Имеют следующий вид:

, .

Начальные условия для заданной системы .

Время слежения .

Задающее воздействие в виде системы ДУ

Начальные условия для воздействия:

.

Введем расширенный вектор состояния и расширенные матрицы

,

,

.

Тогда новое описание системы имеет вид:

С начальными условиями: .

Решением уравнения Риккати будет матрица:

С н. у.

Тогда оптимальное управление, находится по формуле:

Используя скрипт AKOR_slegenie_na_konech_interval_I_podxod, получили следующие результаты:

Рис.41. Графики решения уравнения Риккати.

Рис.42. Графики коэффициентов регулятора обратной и прямой связи.

Рис.43. Графики фазовых координат.

Рис.44. График управления.

Выводы: На данном этапе была решена задача АКОР-слежения. В качестве отслеживаемого воздействия была взята исходная система, но с другими начальными условиями, поэтому графики фазовых координат отличаются от заданных, но только на начальном участке движения.

5.5 Задача АКОР для отслеживания известного задающего воздействия. II подход (линейный сервомеханизм)

Система задана в виде:

Матрицы Заданы в пункте 5.1.1.

Весовые матрицы И Имеют следующий вид:

, .

Начальные условия для заданной системы .

Задающее воздействие имеет вид:

, .

Время слежения

Введем вспомогательную вектор-функцию , ДУ которой определяется

,

,

НУ определяются из соотношения

Зная закон изменения и , можно определить управление:

.

Используя скрипт AKOR_slegenie_na_konech_interval_II_podxod, получили следующие результаты:

Рис.45. Графики решения уравнения Риккати.

Рис.46. График задающего воздействия.

Рис.47. Графики коэффициентов регулятора обратной и прямой связи.

Рис.48. Графики фазовых координат.

Рис.49. График управления.

Выводы: На данном этапе была решена задача построения линейного сервомеханизма. В качестве отслеживаемого воздействия была задана экспоненциальная функция. Анализируя выше приведенные графики, можно сказать, что все состояния заданной системы, особенно первая фазовая координата, отслеживается с заданной точностью.

5.6 Задача АКОР – слежения со скользящими интервалами

Пусть интервал времени является объединением нескольких отрезков. Известно некоторое задающее воздействие заданное аналитическим выражением, причем информация о задающем сигнале на следующем отрезке времени поступает только в конце предыдущего. Таким образом, зная задающий сигнал только на одном отрезке времени, мы будем синтезировать управление на этом отрезке.

Разобьем весь интервал на 3 равных отрезка.

Данная задача похожа на задачу отслеживания известного задающего воздействия, заданного аналитическим выражением, но с некоторыми изменениями:

1. Поскольку в уравнение Риккати относительно матрицы входят только параметры системы и функционала качества, то решать его будем один раз на первом отрезке, так как на остальных отрезках решение будет иметь тот же вид, но будет смещено по времени:

2. Начальными условиями для системы на каждом отрезке будет точка, в которую пришла система на предыдущем отрезке:

3. Вектор необходимо пересчитывать на каждом отрезке.

4. В остальном данная задача аналогична задаче построения линейного сервомеханизма (пункт 5.5).

Используя скрипт AKOR_slegenie_so_skolz_intervalami_Modern, получили следующие результаты:

Рис.50. Графики решения уравнения Риккати.

Рис.51. Графики фазовых координат.

Рис.52. График управления.

Выводы: при сравнении полученных результатов, можно сказать, что различия в фазовых координатах при наличии трех участков и при наличии одного участка несущественные. Если сравнивать скорость вычислений и используемые ресурсы, то скорость увеличивается почти в 3 раза, а памяти требуется в 3 раза меньше для решения поставленной задачи. В точках соединения участков наблюдаются скачки, связанные с тем, что требуется значительные затраты на управление, но для первой координаты этот скачок незначительный.

6. Синтез наблюдателя полного порядка

Наблюдателями называются динамические устройства, которые позволяют по известному входному и выходному сигналу системы управления получить оценку вектора состояния. Причем ошибка восстановления .

Система задана в виде:

Начальные условия для заданной системы .

Матрицы Заданы в пункте 5.1.1.

Весовые матрицы И Имеют следующий вид:

, .

Построим наблюдатель полного порядка и получим значения наблюдаемых координат таких, что:

В качестве начальных условий для наблюдателя выберем нулевые н. у.:

Ранг матрицы наблюдаемости:

– матрица

Наблюдаемости.

.

.

Т. е. система является наблюдаемой.

Коэффициенты регулятора:

,

Тогда

Собственные значения матрицы :

Коэффициенты наблюдателя выберем из условия того, чтобы наблюдатель был устойчивым, и ближайший к началу координат корень матрицы лежал в 3 – 5 раз левее, чем наиболее быстрый корень матрицы . Выберем корни матрицы

Коэффициенты матрицы наблюдателя:

.

Используя скрипт Sintez_nablyud_polnogo_poryadka, получили следующие результаты:

Рис.53. Графики решения уравнения Риккати.

Рис.54. Графики фазовых координат.

Рис.55. Графики управлений.

Выводы : Так как система является полностью наблюдаема и полностью управляема, то спектр матрицы может располагаться произвольно. Перемещая собственные значения матрицы левее, относительно собственных значений матрицы мы улучшаем динамику системы, однако, наблюдатель становится более чувствителен к шумам.

Литература

1. Методы классической и современной теории автоматического управления: Учебник в 5 – и т. Т.4: Теория оптимизации систем автоматического управления / Под ред. Н. Д. Егупова. – М.: Изд-во МГТУ им. Н. Э. Баумана, 2004. – 748 с.

2. Краснощеченко В. И.: Методическое пособие: “Методы теории оптимального управления”.

Приложение.

PlotTimeFrHaract. m

Clc

Clear all

Close all

B1 = 9;

B0 = 5;

A4 = 0.1153;

A3 = 1.78;

A2 = 3.92;

A1 = 14.42;

A0 = 8.583;

% syms s w

% W_s_chislit = b1 * s + b0;

% W_s_znamen = s * (a4 * s^4 + a3 * s^3 + a2 * s^2 + a1 * s + a0);

%

% W_s_obj = W_s_chislit/W_s_znamen;

%A_w = collect(simplify(abs(subs(W_s_obj, s, i*w))))

%———————-Построение АЧХ————————————-%

Figure(‘Name’, ‘[0,10]’);

W = 0 : 0.01 : 10;

A_w = sqrt((b0^2 + b1^2.*w.^2)./((-a1*w.^2+a3*w.^4).^2+(a0*w-a2*w.^3+a4*w.^5).^2));

Plot(w, A_w,’k’, ‘LineWidth’, 2);

Grid on

Xlabel(‘w’)

Ylabel(‘A(w)’)

Title(‘Function ACHX(w)’)

%————————————————————————-%

R_ch = roots([b1 b0])

R_zn = roots([a4 a3 a2 a1 a0 0])

%———————-Построение ФЧХ————————————-%

Figure(‘Name’, ‘[0,100]’);

W = 0 : 0.01 : 100;

Fi_w = (atan(w/0.5556)-atan(w/0)-atan(w/13.5832)-atan((w-2.7677)/0.5850)…

-atan((w+2.7677)/0.5850) – atan(w/(0.6848)))*180/pi;

Plot(w, fi_w, ‘k’, ‘LineWidth’, 2);

Grid on

Xlabel(‘w’)

Ylabel(‘fi(w)’)

Title(‘Function FCHX(w)’)

%————————————————————————-%

%———————-Построение АФЧХ————————————%

Figure(‘Name’, ‘[0,100]’);

W = 0 : 0.01 : 100;

A_w = sqrt((b0^2 + b1^2.*w.^2)./((-a1*w.^2+a3*w.^4).^2+(a0*w-a2*w.^3+a4*w.^5).^2));

Fi_w = (atan(w/0.5556)-atan(w/0)-atan(w/13.5832)-atan((w-2.7677)/0.5850)…

-atan((w+2.7677)/0.5850) – atan(w/(0.6848)));

Polar(fi_w, A_w, ‘k’);

Grid on

Xlabel(‘Re(W(jw))’)

Ylabel(‘Im(W(jw))’)

Title(‘Function AFCHX(fi_w, A_w)’)

%————————————————————————-%

%———————-Построение ЛАЧХ————————————%

Figure(‘Name’, ‘[0,100]’);

W = -100 : 0.01 : 100;

LA_w = 20*log(sqrt((b0^2 + b1^2.*w.^2)./((-a1*w.^2+a3*w.^4).^2+(a0*w-a2*w.^3+a4*w.^5).^2)));

Plot(w, LA_w,’k’, ‘LineWidth’, 2);

Grid on

Xlabel(‘w’)

Ylabel(‘L(w)’)

Title(‘Function L(w)’)

%————————————————————————-%

%———————-Построение ФАЧХ————————————%

%————————————————————————-%

%———————-Построение h(t)————————————%

Figure(‘Name’, ‘[0,50]’);

T = 0 : 0.01 : 50;

H_t = 0.0024 * exp(-13.5832.*t) – 0.2175 * exp(-0.6848.*t)…

+ 0.1452 * exp(-0.5850.*t).* cos(2.7677.*t)…

– 0.2217 * exp(-0.5850.*t).* sin(2.7677.*t)…

+ 0.5825 .* t + 0.0699;

Plot(t, h_t, ‘k’, ‘LineWidth’, 2);

Grid on

Xlabel(‘t’)

Ylabel(‘h(t)’)

Title(‘Function h(t)’)

%————————————————————————-%

%———————-Построение k(t)————————————%

Figure(‘Name’, ‘[0,50]’);

T = 0 : 0.01 : 50;

K_t = – 0.0329 * exp(-13.5832.*t) + 0.1489 * exp(-0.6848.*t)…

– 0.6986 * exp(-0.5850.*t).* cos(2.7677.*t)…

– 0.2721 * exp(-0.5850.*t).* sin(2.7677.*t)…

+ 0.5826;

Plot(t, k_t, ‘k’, ‘LineWidth’, 2);

Grid on

Xlabel(‘t’)

Ylabel(‘k(t)’)

Title(‘Function k(t)’)

%————————————————————————-%

X1=tf([b1 b0],[a4 a3 a2 a1 a0 0]);

Ltiview(x1)

ProstranstvoSostoyanii. m

Clc

Clear all

%format rational

B1 = 9;

B0 = 5;

A5 = 0.1153;

A4 = 1.78;

A3 = 3.92;

A2 = 14.42;

A1 = 8.583;

A0 = 0;

%1. Матрица Фробениуса

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

A=[0 1 0 0 0;

0 0 1 0 0;

0 0 0 1 0;

0 0 0 0 1;

0 – a1/a5 – a2/a5 – a3/a5 – a4/a5]

B=[0; 0; 0; 0; 1/a5]

C=[b0 b1 0 0 0]

%Проверка

Syms s

W_s = collect(simplify(C*(s.*eye(5)-A)^(-1)*B),s)

Pretty(W_s)

%2. Параллельная декомпозиция

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

B1 = b1/a5;

B0 = b0/a5;

S1 = 0;

S2 = -6615/487;

S3 = -1022/1747 + 4016/1451*i;

S4 = -1022/1747 – 4016/1451*i;

S5 = -415/606;

Alfa = real(s3);

Beta = imag(s3);

Syms s A B C D E

W_s_etal = collect(((b1*s+b0)/((s-s1)*(s-s2)*((s+alfa)^2+beta^2)*(s-s5))),s)

%pretty(W_s_etal)

Slag_1 = simplify(collect(A*(s-s2)*((s+alfa)^2+beta^2)*(s-s5),s));

Slag_2 = simplify(collect(B*(s-s1)*((s+alfa)^2+beta^2)*(s-s5),s));

Slag_3 = simplify(collect(C*(s-s1)*((s+alfa)^2+beta^2)*(s-s2),s));

Slag_4 = simplify(collect((D*s+E)*(s-s1)*(s-s2)*(s-s5),s));

Chislit_W_s =collect(Slag_1 + Slag_2 + Slag_3 + Slag_4,s);

%Решение системы линейных уравнений

MS =double( [1 1 1 1 0;

6753029497/515578134 -513659/1058682 10560977/850789 4210795/295122 1;

77456808434995506239663107/126764366837761533378822144 1874500571398143988939141/260296441145300889894912 -3300780600401725219142291/418364246989311991349248 915075/98374 4210795/295122;

26189071674868424275768861465/253528733675523066757644288 2853037197681682345182805/520592882290601779789824 45476725452203201718998205/418364246989311991349248 0 915075/98374;

6290947020888109571128085025/84509577891841022252548096 0 0 0 0])

PCH = [0; 0; 0; b1; b0];

Koeff = MS^(-1)*PCH

%Проверка

MS*[Koeff(1);Koeff(2);Koeff(3);Koeff(4);Koeff(5)];

Slag_1 = simplify(collect(Koeff(1)*(s-s2)*((s+alfa)^2+beta^2)*(s-s5),s));

Slag_2 = simplify(collect(Koeff(2)*(s-s1)*((s+alfa)^2+beta^2)*(s-s5),s));

Slag_3 = simplify(collect(Koeff(3)*(s-s1)*((s+alfa)^2+beta^2)*(s-s2),s));

Slag_4 = simplify(collect((Koeff(4)*s+Koeff(5))*(s-s1)*(s-s2)*(s-s5),s));

Chislit_W_s =collect((Slag_1 + Slag_2 + Slag_3 + Slag_4),s);

Znamena_W_s = collect((s-s1)*(s-s2)*((s+alfa)^2+beta^2)*(s-s5),s);

W_s = collect(simplify(Koeff(1)/(s-s1)+Koeff(2)/(s-s2)+(Koeff(4)*s+Koeff(5))/((s+alfa)^2+beta^2)+Koeff(3)/(s-s5)),s)

Pretty(W_s)

%Расчет матриц состояния

A = [s1 0 0 0 0;

0 s2 0 0 0 ;

0 0 0 1 0;

0 0 -(alfa^2+beta^2) -2*alfa 0;

0 0 0 0 s5]

B = [Koeff(1); Koeff(2); 0; 1; Koeff(3)]

C = [1 1 Koeff(5) Koeff(4) 1]

%Проверка

W_s = collect(simplify(C*(s.*eye(5)-A)^(-1)*B),s)

Pretty(W_s)

%ВСЕ ПОДСЧИТАНО ВЕРНО!!!

SimplexMetod2.m

Function SimplexMetod2

Clc

Clear all

Close all

Format short

%%%%%%%%%%%%%%%%%%%%%%%ВВОДИМЫЕ ДАННЫЕ%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Матрицы системы

A = [0 2;

-3 0];

B = [0; 2];

% Координаты начальной и конечной точки

X_0 = [14; 0];

X_N = [0; 0];

% Ограничение на управление

U_m = -3;

U_p = 5;

Eps = 1e-10;% погрешность сравнения с нулем

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

N = 195;% число шагов

%h = t1/N;% шаг дискретизации

H = 0.0162;

Alfa = 1;

A = 0;

B = 0;

%t1 = 796/245;% время перехода в конечное состояние

N = size(A);

N = n(1);% порядок системы

% Нахождение матричного экспоненциала

Syms s t

MatrEx = ilaplace((s*eye(n)-A)^(-1));

MatrEx_B = MatrEx*B;

% Вычисление матриц F и G

F = subs(MatrEx, t, h);

G = double(int(MatrEx_B, t, 0, h));

%%%%%%%%%%ФОРМИРОВАНИЕ ЗАДАЧИ БЫСТРОДЕЙСТВИЯ КАК ЗАДАЧИ%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%ЛИНЕЙНОГО ПРОГРАММИРОВАНИЯ%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

For index = 1 : 1e+10

% Вычисление правой части

PravChast = X_N – F^N * X_0;

% Вычисление произведения F на G

FG = zeros(n, N);% формирование матрицы для хранения данных

For j = 1 : n

For i = 0 : N – 1

Fg = F^(N-i-1) * G;

If PravChast(j) < 0

Fg = – fg;

End

FG(j, i+1) = fg(j);

End

End

% Построение z-строки

Z_stroka = zeros(1, 4*N+n+2);% формирование матрицы для хранения данных

% Первый элемент z-строки

Z_stroka(1) = 1;

% Суммирование правых частей

For j = 1 : n

Z_stroka(4*N+n+2) = z_stroka(4*N+n+2) + abs(PravChast(j));

End

% Формирование элементов z-строки между 1-м и последним элементами

%при 2N небазисных переменных, т. е. при управлениях

For i = 2 : 2 : 2 * N

For j = 1 : n

Z_stroka(i) = z_stroka(i) + FG(j, i/2);

End

For j = 1 : n

Z_stroka(i+1) = z_stroka(i+1) – FG(j, i/2);

End

End

% Формирование симплекс-таблицы

CT = zeros(n+2*N+1, 4*N+n+2);

% Построение симплекс-таблицы начиная с z-строки

CT(1,:) = z_stroka(1,:);

% Формирование R-строк в симплекс-таблице

For j = 2 : n + 1

% Формирование правой части в R-строках

CT(j, 4*N+n+2) = abs(PravChast(j-1));

% Формирование элементов R-строк между 1-м и последним элементами

%при 2N небазисных переменных, т. е. при управлениях

For i = 2 : 2 : 2 * N

CT(j, i) = FG(j-1, i/2);

CT(j, i+1) = – FG(j-1, i/2);

End

End

% Формирование S-строк в симплекс-таблице

L = 2;

For j = n + 2 : 2 : n + 2 * N + 1

% Формирование правой части в S-строках

CT(j, 4*N+n+2) = u_p;

CT(j+1, 4*N+n+2) = abs(u_m);

% Формирование элементов S-строк между 1-м и последним элементами

%при 2N небазисных переменных, т. е. при управлениях

CT(j, l : l+1) = [1 -1];

CT(j+1, l : l+1) = [-1 1];

L = l + 2;

End

% Формирование базиса в симплекс-таблице, т. е коэффициентов, стоящих при

%базисных переменных от 2N небазисных переменных до правой части (до 4*N+n+1)

CT(2 : n+2*N+1, 2*N+2 : 4*N+n+1) = eye(n+2*N, n+2*N);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%

РЕШЕНИЕ ЗАДАЧИ БЫСТРОДЕЙСТВИЯ%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%СИМПЛЕКС-МЕТОДОМ%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Цикл смены базисных переменных

Nn = size(find(CT(1,2:2*N+1) >= eps));

While nn > 0

[znach, N_stolb] = max(CT(1, 2 : 2*N+1));

N_stolb = N_stolb + 1; % т. к. при небазисн. перемен.

PravChast = CT(:, 4*N+n+2);

For j = 2 : n + 2 * N + 1

if CT(j, N_stolb) > 0

PravChast(j) = PravChast(j) / CT(j, N_stolb);

else

PravChast(j) = inf;

end

End

[znach, N_str] = min(PravChast(2 : n+2*N+1));

N_str = N_str + 1;

% Формирование матрицы перехода B

B = eye(n+2*N+1, n+2*N+1);

B(:, N_str) = CT(:, N_stolb);

% Обращение матрицы B

RE = B(N_str, N_str);

For j = 1 : n + 2 * N + 1

If j == N_str

B(j, N_str) = 1 / RE;

Else

B(j, N_str) = – B(j, N_str) / RE;

End

End

%B = inv(B);

% Получение новой симплекс таблицы

CT = B * CT;

Nn = size(find(CT(1,2:2*N+1) >= eps));

End

U = zeros(1,N);

% Формирование управления

For j = 2 : n + 2 * N + 1

For i = 2 : 2 * N + 1

If CT(j, i) >= eps

If mod(i, 2) < eps

U(i/2) = CT(j, 4*N+n+2);

Else

U((i-1)/2) = – CT(j, 4*N+n+2);

End

End

End

End

% Формирование x1 и x2

X = zeros(n, N);

X(:, 1) = F * X_0 + G * u(1);

For i = 2 : N

X(:, i) = F * X(:, i-1) + G * u(i);

End

% Объединение с начальными условиями

X1 = [X_0(1) X(1, :)];

X2 = [X_0(2) X(2, :)];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% проверка на окончание выбора количества шагов

XX = [X_0 X];

% Вычисление нормы вектора состояния

NormaXX = norm(XX(:,N))

% Вычисление значения переменной R

R = abs(X_N – F^N * X_0) – FG * u’;

R = R’;

Z = sum(R);

% Погрешность приближения к точному решению

Pogresh = 0.3;

If (normaXX < pogresh)

N_opt = N;

Break;

Else

If (z > h)

If a == 1

Alfa = ceil(alfa/2);

End

N = N + alfa;

A = 0;

B = 1;

Else

If b == 1

Alfa = ceil(alfa/2);

End

N = N – alfa;

A = 1;

B = 0;

End

End

T_perevoda = N * h;

End

N_opt

H

T_perevoda

%%%%%%%%%%%%%%%%%%%%ОФОРМЛЕНИЕ ПОЛУЧЕННЫХ РЕЗУЛЬТАТОВ%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%В ГРАФИЧЕСКОМ ВИДЕ%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Построение графика x1(t);

Figure(1)

T = (0 : 1 : length(X1)-1) * h;

Plot(t, X1, ‘b’, ‘LineWidth’, 2);

Hl=legend(‘x_1(t)’);

Set(hl, ‘FontName’, ‘Courier’);

Xlabel(‘t, cek’); ylabel(‘x_1(t)’);

Grid on

% Построение графика x2(t);

Figure(2)

T = (0 : 1 : length(X2)-1) * h;

Plot(t, X2, ‘b’, ‘LineWidth’, 2);

Hl=legend(‘x_2(t)’);

Set(hl, ‘FontName’, ‘Courier’);

Xlabel(‘t, cek’); ylabel(‘x_2(t)’);

Grid on

% Построение графика x2 = x2(x1);

Figure(3)

Plot(X1, X2, ‘m’, ‘LineWidth’, 2);

Hl=legend(‘x_2 = x_2(x_1)’);

Set(hl, ‘FontName’, ‘Courier’);

Xlabel(‘x_1(t)’); ylabel(‘x_2(x_1(t))’);

Grid on

% Построение графика u(t)

Figure(4)

T = (0 : 1 : length(u)-1) * h;

Plot(t, u, ‘r’, ‘LineWidth’, 2);

Hl=legend(‘u(t)’);

Set(hl, ‘FontName’, ‘Courier’);

Xlabel(‘t, cek’); ylabel(‘u(t)’);

Grid on

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Optimal_L_problem_moments. m

Clc

Close all

Clear all

Format long

% ————————————————————————%

B_0 = 5;

B_1 = 9;

% Укороченная система данного объекта

A_5 = 0.1153;

A_4 = 1.78;

A_3 = 3.92;

A_2 = 14.42;

A_1 = 8.583;

A_0 = 0;

% ————————————————————————%

% Приведение системы

B0 = b_0/a_5;

B1 = b_1/a_5;

A5 = a_5/a_5;

A4 = a_4/a_5;

A3 = a_3/a_5;

A2 = a_2/a_5;

A1 = a_1/a_5;

A0 = a_0/a_5;

% ————————————————————————%

% Порядок системы

Poryadok = 5;

% Начальные и конечные условия относительно вектора Y

Y_0 = [3 2 1 5]’;

Y_T = [0 -1 0 3]’;

% Конечное время перехода

T = 3;

% Матрица перехода от Н. У. Y к Н. У. X

B_ = [b0 b1 0 0 0;

0 b0 b1 0 0;

0 0 b0 b1 0;

0 0 0 b0 b1];

% ————————————————————————%

% ————————————————————————%

% Начальные условия для упорядоченной системы

X_0 = B_’ * inv(B_ * B_’) * Y_0

X_T = B_’ * inv(B_ * B_’) * Y_T

% ————————————————————————%

% ————————————————————————%

% Представление системы в пространстве состояний

A = [0 1 0 0 0;

0 0 1 0 0;

0 0 0 1 0

0 0 0 0 1;

-a0 – a1 – a2 – a3 – a4]

B = [0; 0; 0; 0; 1]

C = [b0 b1 0 0 0]

% ————————————————————————%

% ————————————————————————%

% Вычисление матричной экспоненты

Syms s t

MatrEx = simplify (vpa(ilaplace(inv(s*eye(5) – A)), 50))

% ————————————————————————%

RETURN = 1;

While RETURN == 1

Disp(‘L – проблема моментов в пространстве вход-выход: 1’)

Disp(‘L – проблема моментов в пространстве состояний : 2’)

Reply = input(‘Выберете метод решения [1 или 2]: ‘, ‘s’);

Switch reply

case ‘1’

disp(‘L – проблема моментов в пространстве вход-выход’)

% ————————L – проблема моментов—————————%

% ———————-в пространстве вход-выход————————-%

% ————————————————————————%

% Передаточная функция

W_obj_s = 1/(a5*s^5 + a4*s^4 + a3*s^3 + a2*s^2 + a1*s + a0);

% Полюса передаточной функции

Polyusa_TF = roots([a5 a4 a3 a2 a1 a0]);

% ИПФ

K_t = simplify (vpa (ilaplace(1 / (a5*s^5 + a4*s^4 + a3*s^3 + …

a2*s^2 + a1*s + a0)),50))

% K_t = vpa(K_t,6)

%... ————————————————————————%

% Составление матрицы Вронского

For i = 1 : poryadok

Matrix_Vron (i, 1) = diff (exp (polyusa_TF(1) *t), t, i – 1);

Matrix_Vron (i, 2) = diff (exp (polyusa_TF(2) *t), t, i – 1);

Matrix_Vron (i, 3) = diff (exp (real(polyusa_TF(3))*t) * …

Cos(imag(polyusa_TF(3))*t), t, i – 1);

Matrix_Vron (i, 4) = diff (exp (real(polyusa_TF(4))*t) * …

Sin(imag(polyusa_TF(4))*t), t, i – 1);

Matrix_Vron (i, 5) = diff (exp (polyusa_TF(5) *t), t, i – 1);

End

% Матрица Вронского при t = 0;

Matrix_Vron_t_0 = double(subs(Matrix_Vron, t,0));

% Матрица Вронского при t = T;

T = 3;

Matrix_Vron_t_T = double(subs(Matrix_Vron, t,T));

% vpa(Matrix_Vron_t_0,6)

% ————————————————————————%

% ————————————————————————%

% Определение неизвестных коэффициентов C

C_ = inv(Matrix_Vron_t_0) * X_0;

% ————————————————————————%

% ————————————————————————%

% Нахождение моментных функций

K_Tt_1 = subs (K_t, t, T – t);

K_Tt = diff (K_t);

K_Tt_2 = subs (K_Tt, t, T – t);

K_Ttt = diff (K_Tt);

K_Tt_3 = subs (K_Ttt, t, T – t);

K_Tttt = diff (K_Ttt);

K_Tt_4 = subs (K_Tttt, t, T – t);

K_Ttttt = diff (K_Tttt);

K_Tt_5 = subs (K_Ttttt, t, T – t);

H1_Tt = K_Tt_1

H2_Tt = K_Tt_2

H3_Tt = K_Tt_3

H4_Tt = K_Tt_4

H5_Tt = K_Tt_5

% ————————————————————————%

% ————————————————————————%

% Нахождение моментов

For i = 1 : poryadok

Matrix_a(i) = X_T(i) – C_’ * Matrix_Vron_t_T(i,:)’;

End

Matrix_a = Matrix_a’

% ————————————————————————%

% ————————————————————————%

% ————————————————————————%

RETURN = 2;

case ‘2’

disp(‘L – проблема моментов в пространстве состояний’)

% ————————L – проблема моментов—————————%

% ———————-в пространстве состояний————————–%

% ————————————————————————%

Matr_Ex_T = subs(MatrEx, t, T);

% ————————————————————————%

% ————————————————————————%

% Нахождение моментов

For i = 1 : poryadok

Matrix_a(i) = X_T(i) – Matr_Ex_T(i,:) * X_0;

End

Matrix_a = Matrix_a’

% ————————————————————————%

% ————————————————————————%

% Нахождение моментных функций

Matr_Ex_Tt = subs(MatrEx, t, T – t);

H_Tt = vpa(expand(simplify(Matr_Ex_Tt * B)),50);

H1_Tt = h_Tt(1)

H2_Tt = h_Tt(2)

H3_Tt = h_Tt(3)

H4_Tt = h_Tt(4)

H5_Tt = h_Tt(5)

% ————————————————————————%

% ————————————————————————%

% ————————————————————————%

RETURN = 2;

Otherwise

disp(‘Неизвестный метод.’)

RETURN = 1;

End

End

% h1_Tt = vpa(h1_Tt,6)

% h2_Tt = vpa(h2_Tt,6)

% h3_Tt = vpa(h3_Tt,6)

% h4_Tt = vpa(h4_Tt,6)

% h5_Tt = vpa(h5_Tt,6)

% ————————————————————————%

% ——–Нахождение управления и вычисление минимальной энергии———-%

% ————————————————————————%

Syms ks1 ks2 ks3 ks4 ks5

% ————————————————————————%

% Формирование функционала

D_v_2 = vpa (simplify (int ((ks1*h1_Tt + ks2*h2_Tt + ks3*h3_Tt + …

ks4*h4_Tt + ks5*h5_Tt)^2, t, 0, T)), 50);

% Выражаем ks1 через остальные

Ks1 = vpa ((1 – ks2*Matrix_a(2) – ks3*Matrix_a(3) – …

ks4*Matrix_a(4) – ks5*Matrix_a(5))/Matrix_a(1), 50);

% Подставляем в функционал ks1

D_v_2 = vpa (expand (subs (d_v_2, ks1)), 50);

% Находим частные производные по ksi

Eq_1= diff(d_v_2, ks2);

Eq_2= diff(d_v_2, ks3);

Eq_3= diff(d_v_2, ks4);

Eq_4= diff(d_v_2, ks5);

% Решаем СЛАУ относительно ksi

Ksi= solve(eq_1, eq_2, eq_3, eq_4);

% Полученные значения ksi

Ks2= double(ksi. ks2)

Ks3= double(ksi. ks3)

Ks4= double(ksi. ks4)

Ks5= double(ksi. ks5)

Ks1 = double(vpa ((1 – ks2*Matrix_a(2) – ks3*Matrix_a(3) – ks4*Matrix_a(4) – …

Ks5*Matrix_a(5))/Matrix_a(1), 50))

% ————————————————————————%

% ————————————————————————%

% Проверка условия полученного результата

ks1*Matrix_a(1) + ks2*Matrix_a(2) + ks3*Matrix_a(3) + …

ks4*Matrix_a(4) + ks5*Matrix_a(5)

% ————————————————————————%

% ————————————————————————%

% Вычисление управления и минимальной энергии

D_v_2 = vpa (simplify (int ((ks1*h1_Tt + ks2*h2_Tt + ks3*h3_Tt + …

ks4*h4_Tt + ks5*h5_Tt)^2, t, 0, T)), 50)

% d_v_2 = double(d_v_2)

Gamma_v_2 = 1/d_v_2

% gamma_v_2 = double(gamma_v_2)

U = vpa (expand(simplify(gamma_v_2 * (ks1*h1_Tt + ks2*h2_Tt + ks3*h3_Tt + …

ks4*h4_Tt + ks5*h5_Tt))), 50)

% u = vpa(u,6)

U_0 = subs(u, t,0)

U_T = subs(u, t,T)

Ezplot(u, [0 T], 1)

Hl=legend(‘u(t)’);

Set(hl, ‘FontName’, ‘Courier’);

Title (‘u(t)’);

Xlabel(‘t’)

Grid on

% ————————————————————————%

% ————————————————————————%

% Нахождения X

% Вычисление матричной экспоненты

MatrEx = simplify (vpa(ilaplace(inv(s*eye(5) – A)), 50));

Syms t tay

X_svob = MatrEx * X_0;

X_vinyg = int ((subs(MatrEx, t, t – tay))*B*(subs (u, t, tay)), tay, 0,t);

X_real = X_svob + X_vinyg;

Save Sostoyaniya X_real u

X_real = vpa (expand (simplify(X_real)), 50)

X_real_0 = double(subs (X_real, t, 0))

X_real_T = double(subs (X_real, t, T))

% Погрешность X

Delta_X_T = double(vpa(X_T – X_real_T, 50))

Delta_X_0 = double(vpa(X_0 – X_real_0, 50))

% Нахождение Y

For i = 1 : poryadok – 1

Y_real(i) = B_(i,:) * X_real;

End

Y_real = vpa (expand(simplify(Y_real’)), 50)

Y_real_0 = double(subs (Y_real, t, 0))

Y_real_T = double(subs (Y_real, t, T))

% Погрешность Y

Delta_Y_T = double(vpa(Y_T – Y_real_T, 50))

Delta_Y_0 = double(vpa(Y_0 – Y_real_0, 50))

% ————————————————————————%

% ————————————————————————%

% Вычисление max значений для задачи АКОР

H = 0.01;

Tic

Tt = 0 : h : T;

For i = 1 : poryadok

X_max(i) = max(abs(subs(X_real(i),t, tt)));

End

U_max = max(abs(subs(u, t,tt)));

Toc

Save Sostoyaniya X_max U_max

% ————————————————————————%

% ————————————————————————%

% Построение результатов X(t)

Ezplot (X_real(1), [0 T],2)

Title (‘x_1(t)’);

Grid on

Ezplot (X_real(2), [0 T],3)

Title (‘x_2(t)’);

Grid on

Ezplot (X_real(3), [0 T],4)

Title (‘x_3(t)’);

Grid on

Ezplot (X_real(4), [0 T],5)

Title (‘x_4(t)’);

Grid on

Ezplot (X_real(5), [0 T],6)

Title (‘x_5(t)’);

Grid on

% Построение результатов Y(t)

Ezplot (Y_real(1), [0 T],7)

Title (‘y_1(t)’);

Grid on

Ezplot (Y_real(2), [0 T],8)

Title (‘y_2(t)’);

Grid on

Ezplot (Y_real(3), [0 T],9)

Title (‘y_3(t)’);

Grid on

Ezplot (Y_real(4), [0 T],10)

Title (‘y_4(t)’);

Grid on

% ————————————————————————%

Gramian_Uprav. m

Clc

Close all

Clear all

Format long

% ————————————————————————%

B_0 = 5;

B_1 = 9;

% Укороченная система данного объекта

A_5 = 0.1153;

A_4 = 1.78;

A_3 = 3.92;

A_2 = 14.42;

A_1 = 8.583;

A_0 = 0;

% ————————————————————————%

% Приведение системы

B0 = b_0/a_5;

B1 = b_1/a_5;

A5 = a_5/a_5;

A4 = a_4/a_5;

A3 = a_3/a_5;

A2 = a_2/a_5;

A1 = a_1/a_5;

A0 = a_0/a_5;

% ————————————————————————%

% Порядок системы

Poryadok = 5;

% Начальные и конечные условия относительно вектора Y

Y_0 = [3 2 1 5]’;

Y_T = [0 -1 0 3]’;

% Конечное время перехода

T = 3;

% Матрица перехода от Н. У. Y к Н. У. X

B_ = [b0 b1 0 0 0;

0 b0 b1 0 0;

0 0 b0 b1 0;

0 0 0 b0 b1];

% ————————————————————————%

% ————————————————————————%

% Начальные условия для упорядоченной системы

X_0 = B_’ * inv(B_ * B_’) * Y_0

X_T = B_’ * inv(B_ * B_’) * Y_T

% ————————————————————————%

% ————————————————————————%

% Представление системы в пространстве состояний

A = [0 1 0 0 0;

0 0 1 0 0;

0 0 0 1 0

0 0 0 0 1;

-a0 – a1 – a2 – a3 – a4];

B = [0; 0; 0; 0; 1];

C = [b0 b1 0 0 0];

% ————————————————————————%

% ————————————————————————%

% Вычисление матричной экспоненты

Syms s t

MatrEx = simplify (vpa(ilaplace(inv(s*eye(5) – A)), 50));

MatrEx_T = vpa(subs(MatrEx, t, T),50);

MatrEx_Tt = vpa(subs(MatrEx, t, T-t),50);

% ————————————————————————%

% ————————————————————————%

% Вычисление матрицы управляемости

M_c = [B A*B A^2*B A^3*B A^4*B]

Rank_M_c = rank(M_c); %ранк = 5 – система управляема

% ————————————————————————%

% ————————————————————————%

% Вычисление грамиана управляемости

W_Tt = double(vpa(simplify(int(MatrEx_Tt*B*B’*MatrEx_Tt’,t,0,T)),50))

% ————————————————————————%

% ————————————————————————%

% Формирование управления

U = vpa(expand(simplify(B’*MatrEx_Tt’*inv(W_Tt)*(X_T-MatrEx_T*X_0))),50)

U_0 = subs(u, t,0)

U_T = subs(u, t,T)

U = vpa(u,6)

% ————————————————————————%

Ezplot(u, [0 T], 1)

Title (‘u(t)’);

Xlabel(‘t’)

Grid on

Tt = 0 : 0.01 : T;

U2 = -20.605579750692850622177761310569*exp(-40.749492463732569440253455897187+13.583164154577523146751151965729*t)+19.011167813350479567880663060491*exp(-2.0544534472800777280645828326668+.68481781576002590935486094422228*t)+1.3356706538317879679656856470126*exp(-1.7550088311372150108106250409710+.58500294371240500360354168032368*t)*cos(-8.3032397968812277095785721047505+2.7677465989604092365261907015835*t)+7.2830359327562658520685140088852*exp(-1.7550088311372150108106250409710+.58500294371240500360354168032368*t)*sin(-8.3032397968812277095785721047505+2.7677465989604092365261907015835*t)-8.6096491449877801097840179781687;

U1 = subs(u2, t, tt);

U2 = subs(u, t, tt);

Figure(2)

Plot(tt, u1,’r’,tt, u2,’b’,’LineWidth’,2)

Hl=legend(‘u(t) при решении оптимальной L-проблемы моментов’,’u(t) с использованием грамиана управляемости’);

Set(hl, ‘FontName’, ‘Courier’);

Xlabel(‘t, cek’); ylabel(‘u(t)’);

Title(‘u(t)’)

Grid on

AKOR_stabilizaciya_na_polybeskon_interval. m

Clc

Clear all

Close all

Poryadok = 5;

% ————————————————————————%

B_0 = 5;

B_1 = 9;

% Укороченная система данного объекта

A_5 = 0.1153;

A_4 = 1.78;

A_3 = 3.92;

A_2 = 14.42;

A_1 = 8.583;

A_0 = 0;

% ————————————————————————%

% Приведение системы

B0 = b_0/a_5;

B1 = b_1/a_5;

A5 = a_5/a_5;

A4 = a_4/a_5;

A3 = a_3/a_5;

A2 = a_2/a_5;

A1 = a_1/a_5;

A0 = a_0/a_5;

% ————————————————————————%

% ————————————————————————%

% Представление системы в пространстве состояний

A = [0 1 0 0 0;

0 0 1 0 0;

0 0 0 1 0;

0 0 0 0 1;

-a0 – a1 – a2 – a3 – a4]

B = [0; 0; 0; 0; 1]

C = [b0 b1 0 0 0]

% Начальные условия

X_0 = [10; 0; 6; 4; 8]

%T = 1;

Time = 1;

% ————————————————————————%

% ————————————————————————%

% Получение max значений из файла

Load Sostoyaniya X_max U_max

% ————————————————————————%

% Нахождение элементов матриц Q и R

R(1) = 0.1;

Q(1) = 1/poryadok * r(1) * (U_max)^2 / (X_max(1))^2;

For i = 2 : poryadok

Q(i) = q(1) * (X_max(1))^2 / (X_max(i))^2;

End

Q = diag(q)

R = diag(r)

% Для изменения коэффициентов

% Q(1,1) = Q(1,1);

% Q(2,2) = Q(2,2);

% Q(3,3) = Q(3,3);

% Q(4,4) = Q(4,4);

% Q(5,5) = Q(5,5);

Q(1,1) = Q(1,1)*1e+12;

Q(2,2) = Q(2,2)*1e+8;

Q(3,3) = Q(3,3)*1e+7;

Q(4,4) = Q(4,4)*1e+0;

Q(5,5) = Q(5,5)*1e+2;

R(1,1) = R(1,1);

% ————————————————————————%

% ————————————————————————%

% Решение уравнения Риккати методом диагонализации

P1 = Solve_Riccati_Method_Diag(A, B,Q, R)

% ————————————————————————%

P_nach = zeros(poryadok, poryadok);%+ones(poryadok, poryadok);

% ————————————————————————%

% Решение уравнения Риккати методом обратного интегрирования

P2 = Solve_Riccati_Method_Revers_Integr(A, B,Q, R,Time, poryadok, P_nach)

% ————————————————————————%

% Сравнение расхождения методов

Delta_P = abs(P1-P2)

% Построение графика коэффициентов регулятора

Load Solve_Riccati_Method_Revers_Integr Time_R P N_str

PP = P;

For i = 1 : N_str

P = reshape(PP(i, :), poryadok, poryadok);

K(i, 🙂 = – inv(R)*B’*P;

End

Figure(2)

Plot(Time_R, K(:,1),’-‘,Time_R, K(:,2),’-‘,Time_R, K(:,3),’-‘,Time_R, K(:,4),’-‘,Time_R, K(:,5),’-‘, ‘LineWidth’, 2);

Xlabel(‘t’)

Tit1 = title(‘Коэффициенты обратной связи в прямом времени’);

Set(tit1,’FontName’,’Courier’);

Hl=legend(‘k_1_о_с’,’k_2_о_с’,’k_3_о_с’,’k_4_о_с’,’k_5_о_с’,0);

Set(hl,’FontName’,’Courier’);

Grid on;

% ————————————————————————%

% Решение уравнения Риккати с помощью встроенной функции

% P = vpa(care(A, B,Q, R), 10)

% ————————————————————————%

% ————————————————————————%

% Нахождение коэффициентов регулятора

Disp(‘Коэффициенты регулятора:’)

K1 = – inv(R) * B’ * P1

K2 = – inv(R) * B’ * P2

% ————————————————————————%

% ————————————————————————%

A1_ = A + B * K1;

A2_ = A + B * K2;

% Вычисление матричной экспоненты

Syms s t

MatrEx1 = simplify (vpa(ilaplace(inv(s*eye(5) – A1_)), 50));

MatrEx2 = simplify (vpa(ilaplace(inv(s*eye(5) – A2_)), 50));

% Нахождение координат состояния

X1 = vpa(simplify(MatrEx1 * X_0), 50);

X2 = vpa(simplify(MatrEx2 * X_0), 50);

% Нахождение управления

U1 = vpa(simplify(K1 * X1),50)

U2 = vpa(simplify(K2 * X2),50)

% ————————————————————————%

% Построение u(t) и X(t)

T_sravneniya = 0.2;

Figure(3);

Tt = 0 : 0.01 : T_sravneniya;

Uu1 = subs(u1,t, tt);

Uu2 = subs(u2,t, tt);

Plot(tt, uu1, tt, uu2, ‘LineWidth’, 2)

Title (‘u(t)’);

Xlabel(‘t’)

Hl=legend(‘u(t) – управление’,0);

Set(hl,’FontName’,’Courier’);

Grid on

Ezplot(X1(1), [0 Time], 4)

Hold on

Title (‘x_1(t)’);

Xlabel(‘t’)

Grid on

Ezplot(X1(2), [0 Time], 5)

Title (‘x_2(t)’);

Xlabel(‘t’)

Grid on

Ezplot(X1(3), [0 Time], 6)

Title (‘x_3(t)’);

Xlabel(‘t’)

Grid on

Ezplot(X1(4), [0 Time], 7)

Title (‘x_4(t)’);

Xlabel(‘t’)

Grid on

Ezplot(X1(5), [0 Time], 8)

Title (‘x_5(t)’);

Xlabel(‘t’)

Grid on

Tt = 0 : 0.01 : T_sravneniya;

X21 = subs(X1(1), t, tt);

X22= subs(X1(2), t, tt);

X23= subs(X1(3), t, tt);

X24= subs(X1(4), t, tt);

X25= subs(X1(5), t, tt);

Save Sravnenie_stabilizacii_1 X21 X22 X23 X24 X25 uu1

AKOR_stabilizaciya_na_konech_interval. m

Clc

Clear all

Close all

Poryadok = 5;

% ————————————————————————%

B_0 = 5;

B_1 = 9;

% Укороченная система данного объекта

A_5 = 0.1153;

A_4 = 1.78;

A_3 = 3.92;

A_2 = 14.42;

A_1 = 8.583;

A_0 = 0;

% ————————————————————————%

% Приведение системы

B0 = b_0/a_5;

B1 = b_1/a_5;

A5 = a_5/a_5;

A4 = a_4/a_5;

A3 = a_3/a_5;

A2 = a_2/a_5;

A1 = a_1/a_5;

A0 = a_0/a_5;

% ————————————————————————%

% ————————————————————————%

% Представление системы в пространстве состояний

A = [0 1 0 0 0;

0 0 1 0 0;

0 0 0 1 0

0 0 0 0 1;

-a0 – a1 – a2 – a3 – a4];

B = [0; 0; 0; 0; 1];

C = [b0 b1 0 0 0];

% Начальные условия

X_0 = [10; 0; 6; 4; 8];

Time = 0.2;

% ————————————————————————%

% ————————————————————————%

% Получение max значений из файла

Load Sostoyaniya X_max U_max

% ————————————————————————%

% Нахождение элементов матриц Q и R

% r(1) = 100;

R(1) = 0.1;

Q(1) = 1/poryadok * r(1) * (U_max)^2 / (X_max(1))^2;

For i = 2 : poryadok

Q(i) = q(1) * (X_max(1))^2 / (X_max(i))^2;

End

Q = diag(q);

R = diag(r);

% Для изменения коэффициентов

Q(1,1) = Q(1,1)*1e+12;

Q(2,2) = Q(2,2)*1e+8;

Q(3,3) = Q(3,3)*1e+7;

Q(4,4) = Q(4,4)*1e+0;

Q(5,5) = Q(5,5)*1e+2;

R(1,1) = R(1,1);

% P_prib = eye(poryadok, poryadok);

% P_prib(1,1) = 100;

% P_prib(2,2) = 10;

% % P_prib(3,3) = 1000;

% % P_prib(4,4) = 10;

% % P_prib(5,5) = 1;

% ————————————————————————%

P_nach = zeros(poryadok, poryadok);% + P_prib;

% ————————————————————————%

% Решение уравнения Риккати методом обратного интегрирования

P = Solve_Riccati_Method_Revers_Integr(A, B,Q, R,Time, poryadok, P_nach)

% ————————————————————————%

% Нахождение переменных коэффициентов регулятора

Load Solve_Riccati_Method_Revers_Integr Time_R P N_str

PP = P;

For i = 1 : N_str

P = reshape(PP(i, :), poryadok, poryadok);

K(i, 🙂 = – inv(R)*B’*P;

End

% ————————————————————————%

% Формирование вектора коэффициентов регулятора

% и решения уравнения Риккати в прямом порядке

Load Solve_Riccati_Method_Revers_Integr P

Size(K)

I = 1;

Len_K = length(K(:,1))

For j = len_K : -1 : 1

K_pr(i,:) = K(j,:);

I = i + 1;

End

% ————————————————————————%

% Построение графика переменных коэффициентов регулятора в прямом времени

Figure(2)

Plot(Time_R, K(:,1),’-‘,Time_R, K(:,2),’-‘,Time_R, K(:,3),’-‘,…

Time_R, K(:,4),’-‘,Time_R, K(:,5),’-‘, ‘LineWidth’, 2);

Grid on;

Title(‘K(t)’)

Xlabel(‘t’)

Legend(‘k_1′,’k_2′,’k_3′,’k_4′,’k_5’);

% ————————————————————————%

% ————————————————————————%

For k = 1 : len_K

A_(:,:,k) = A + B * K(k,:);

End

Size(A_);

% ————————————————————————%

% Нахождение фазовых координат

X(:,1) = X_0;

H = 0.01;

Time_X(1) = 0;

For k = 1 : len_K

X(:, k+1) = X(:, k) + h * A_(:,:,k) * X(:, k);

Time_X(k+1) = time_X(k) + h;

End

X(:, k+1) = [];

Time_X(k+1) = [];

% ————————————————————————%

% Нахождение управления

For k = 1 : len_K

U(k) = K_pr(k,:) * X(:,k);

End

% ————————————————————————%

% Построение u(t) и X(t)

Figure(3);

Plot(time_X, u, ‘r-‘, ‘LineWidth’, 2)

Title (‘u(t)’);

Xlabel(‘t’)

Grid on

Figure(4);

Plot(time_X, X(1,:), ‘LineWidth’, 2)

Hold on

Title (‘x_1(t)’);

Xlabel(‘t’)

Grid on

Figure(5);

Plot(time_X, X(2,:), ‘LineWidth’, 2)

Title (‘x_2(t)’);

Xlabel(‘t’)

Grid on

Figure(6);

Plot(time_X, X(3,:), ‘LineWidth’, 2)

Title (‘x_3(t)’);

Xlabel(‘t’)

Grid on

Figure(7);

Plot(time_X, X(4,:), ‘LineWidth’, 2)

Title (‘x_4(t)’);

Xlabel(‘t’)

Grid on

Figure(8);

Plot(time_X, X(5,:), ‘LineWidth’, 2)

Title (‘x_5(t)’);

Xlabel(‘t’)

Grid on

Save Sravnenie_stabilizacii_2 time_X X u

Sravnenie_stabilizacii. m

Close all

Load Sravnenie_stabilizacii_1 X21 X22 X23 X24 X25 uu1

Load Sravnenie_stabilizacii_2 time_X X u

Figure(31);

Plot(time_X, u, time_X, uu1, ‘LineWidth’, 2)

Title (‘u(t)’);

Xlabel(‘t’)

Hl=legend(‘u(t) – управление с перемен. коеф.’,’u(t) – управление с пост. коеф.’);

Set(hl,’FontName’,’Courier’);

Grid on

Figure(41);

Plot(time_X, X(1,:), time_X, X21, ‘LineWidth’, 2)

Hold on

Title (‘x_1(t)’);

Xlabel(‘t’)

Hl=legend(‘x_1(t) – с перемен. коеф.’,’x_1(t) – с пост. коеф.’);

Set(hl,’FontName’,’Courier’);

Grid on

Figure(51);

Plot(time_X, X(2,:), time_X, X22,’LineWidth’, 2)

Title (‘x_2(t)’);

Xlabel(‘t’)

Hl=legend(‘x_2(t) – с перемен. коеф.’,’x_2(t) – с пост. коеф.’);

Set(hl,’FontName’,’Courier’);

Grid on

Figure(61);

Plot(time_X, X(3,:), time_X, X23,’LineWidth’, 2)

Title (‘x_3(t)’);

Xlabel(‘t’)

Hl=legend(‘x_3(t) – с перемен. коеф.’,’x_3(t) – с пост. коеф.’);

Set(hl,’FontName’,’Courier’);

Grid on

Figure(71);

Plot(time_X, X(4,:), time_X, X24,’LineWidth’, 2)

Title (‘x_4(t)’);

Xlabel(‘t’)

Hl=legend(‘x_4(t) – с перемен. коеф.’,’x_4(t) – с пост. коеф.’);

Set(hl,’FontName’,’Courier’);

Grid on

Figure(81);

Plot(time_X, X(5,:), time_X, X25,’LineWidth’, 2)

Title (‘x_5(t)’);

Xlabel(‘t’)

Hl=legend(‘x_5(t) – с перемен. коеф.’,’x_5(t) – с пост. коеф.’);

Set(hl,’FontName’,’Courier’);

Grid on

AKOR_stabilizaciya_pri_vozmusheniyah. m

Clc

Clear all

Close all

Warning off

Poryadok = 5;

% ————————————————————————%

B_0 = 5;

B_1 = 9;

% Укороченная система данного объекта

A_5 = 0.1153;

A_4 = 1.78;

A_3 = 3.92;

A_2 = 14.42;

A_1 = 8.583;

A_0 = 0;

% ————————————————————————%

% Приведение системы

B0 = b_0/a_5;

B1 = b_1/a_5;

A5 = a_5/a_5;

A4 = a_4/a_5;

A3 = a_3/a_5;

A2 = a_2/a_5;

A1 = a_1/a_5;

A0 = a_0/a_5;

% ————————————————————————%

% ————————————————————————%

% Представление системы в пространстве состояний

A = [0 1 0 0 0;

0 0 1 0 0;

0 0 0 1 0

0 0 0 0 1;

-a0 – a1 – a2 – a3 – a4];

B = [0; 0; 0; 0; 1];

C = [b0 b1 0 0 0];

% Начальные условия

X_0 = [10; 0; 6; 4; 8];

Time = 1;

H = 0.01;

% ————————————————————————%

Tic

% ————————————————————————%

% Получение max значений из файла

Load Sostoyaniya X_max U_max

% ————————————————————————%

% Нахождение элементов матриц Q и R

R(1) = 100;

Q(1) = 1/poryadok * r(1) * (U_max)^2 / (X_max(1))^2;

For i = 2 : poryadok

Q(i) = q(1) * (X_max(1))^2 / (X_max(i))^2;

End

Q = diag(q);

R = diag(r);

% Для изменения коэффициентов

Q(1,1) = Q(1,1)*1e+12;

Q(2,2) = Q(2,2)*1e+8;

Q(3,3) = Q(3,3)*1e+7;

Q(4,4) = Q(4,4)*1e+0;

Q(5,5) = Q(5,5)*1e+2;

R(1,1) = R(1,1);

% P_0 = ones(poryadok, poryadok);

% P_0(1,1) = P_0(1,1)*1e12;

% P_0(2,2) = P_0(2,2)*1e8;

% P_0(3,3) = P_0(3,3)*1e7;

% P_0(4,4) = P_0(4,4)*1e0;

% P_0(5,5) = P_0(5,5)*1e2;

% ————————————————————————%

P_nach = zeros(poryadok, poryadok);%+P_0;

% ————————————————————————%

% Решение уравнения Риккати методом обратного интегрирования

P = Solve_Riccati_Method_Revers_Integr(A, B,Q, R,Time, poryadok, P_nach);

Load Solve_Riccati_Method_Revers_Integr_for_slegenie Time_R P N_str

PP = P;

For k = 1 : N_str

P1 = reshape(PP(k, :), poryadok, poryadok);

For i = 1 : poryadok

For j = 1 : poryadok

P2(i, j,k) = P1(i, j);

End

End

End

Size_P = size(P2);

% ————————————————————————%

Tic

% ————————————————————————%

% Получение дискретных значений задающего воздействия в обратном времени

% для нахождения вспомогательной функции q(t)

Vozmyshyayushee_Vozdeistvie_Discrete_Revers(h, 0, Time);

% ————————————————————————%

Load Vozmyshyayushee_Vozdeistvie_Discrete_Revers w_discrete_rev

% ————————————————————————%

Size(w_discrete_rev);

% Начальное значение q(t)

Q = zeros(poryadok,1);

% Интегрирование q(t) в обратном времени

For k = 1 : N_str

Q(:, k+1) = q(:, k) – h * ((P2(:,:,k)*B*inv(R)*B’-A’) * q(:, k) – P2(:,:,k)*w_discrete_rev(:,k));

End

Q(:, k+1) = [];

Size_q = size(q);

% ————————————————————————%

% Нахождение переменных коэффициентов регулятора

For k = 1 : N_str

K_o(k, 🙂 = – inv(R) * B’ * P2(:,:,k);

K_pr(k, 🙂 = – inv(R) * B’;

End

% Формирование вектора коэффициентов регулятора, значений задающего

% воздействия, значений вспомогательной функции в прямом порядке

Size(K_o);

Size(K_pr);

K_pr_p = K_pr;

I = 1;

Len_K = length(K_o(:,1));

For j = len_K : -1 : 1

K_o_p(i,:) = K_o(j,:);

W_discrete(:,i) = w_discrete_rev(:,j);

Q_pr(:, i) = q(:, j);

I = i + 1;

End

% ————————————————————————%

% Построение графика переменных коэффициентов регулятора обратной связи

% в прямом времени

Toc

Figure(3)

Plot(Time_R, K_o(:,1),’-‘,Time_R, K_o(:,2),’-‘,Time_R, K_o(:,3),’-‘,…

Time_R, K_o(:,4),’-‘,Time_R, K_o(:,5),’-‘, ‘LineWidth’, 2);

Xlabel(‘t’)

Tit1 = title(‘Коэффициенты обратной связи в прямом времени’);

Set(tit1,’FontName’,’Courier’);

Hl=legend(‘k_1_о_с’,’k_2_о_с’,’k_3_о_с’,’k_4_о_с’,’k_5_о_с’,0);

Set(hl,’FontName’,’Courier’);

Grid on;

% ————————————————————————%

% ————————————————————————%

% Построение графика переменных коэффициентов регулятора прямой связи

% в прямом времени

Figure(4)

Plot(Time_R, K_pr(:,1),’-‘,Time_R, K_pr(:,2),’-‘,Time_R, K_pr(:,3),’-‘,…

Time_R, K_pr(:,4),’-‘,Time_R, K_pr(:,5),’-‘, ‘LineWidth’, 2);

Xlabel(‘t’)

Tit1 = title(‘Коэффициенты прямой связи в прямом времени’);

Set(tit1,’FontName’,’Courier’);

Hl=legend(‘k_1_п_с’,’k_2_п_с’,’k_3_п_с’,’k_4_п_с’,’k_5_п_с’,0);

Set(hl,’FontName’,’Courier’);

Grid on;

% ————————————————————————%

Tic

% ————————————————————————%

For k = 1 : len_K

A_(:,:,k) = A + B * K_o_p(k,:);

End

Size_A_ = size(A_);

% ————————————————————————%

% ————————————————————————%

% Нахождение фазовых координат

X(:,1) = X_0;

Time_X(1) = 0;

For k = 1 : len_K

X(:, k+1) = X(:, k) + h * (A_(:,:,k) * X(:, k) + B * K_pr_p(k,:) * q_pr(:,k) + w_discrete(:,k));

Time_X(k+1) = time_X(k) + h;

End

X(:, k+1) = [];

Time_X(k+1) = [];

Size_X = size(X);

% ————————————————————————%

% ————————————————————————%

% Нахождение управления

For k = 1 : len_K

U(k) = K_o_p(k,:) * X(:,k) + K_pr_p(k,:) * q_pr(:,k);

End

Size_u = size(u);

% ————————————————————————%

Toc

% Построение u(t) и X(t)

Figure(5);

Plot(time_X, u, ‘r-‘, ‘LineWidth’, 2)

Title (‘u(t)’);

Xlabel(‘t’)

Hl=legend(‘u(t) – управление’,0);

Set(hl,’FontName’,’Courier’);

Grid on

Figure(6);

Plot(time_X, X(1,:),’r-‘, time_X, w_discrete(1,:), ‘LineWidth’, 2)

Hold on

Title (‘x_1(t)’);

Xlabel(‘t’);

Hl=legend(‘X(t) – реальный сигнал’,’w(t) – возмущающее воздействие’,0);

Set(hl,’FontName’,’Courier’);

Grid on

Figure(7);

Plot(time_X, X(2,:),’r-‘, time_X, w_discrete(2,:), ‘LineWidth’, 2)

Title (‘x_2(t)’);

Xlabel(‘t’);

Hl=legend(‘X(t) – реальный сигнал’,’w(t) – возмущающее воздействие’,0);

Set(hl,’FontName’,’Courier’);

Grid on

Figure(8);

Plot(time_X, X(3,:),’r-‘, time_X, w_discrete(3,:), ‘LineWidth’, 2)

Title (‘x_3(t)’);

Xlabel(‘t’);

Hl=legend(‘X(t) – реальный сигнал’,’w(t) – возмущающее воздействие’,0);

Set(hl,’FontName’,’Courier’);

Grid on

Figure(9);

Plot(time_X, X(4,:),’r-‘, time_X, w_discrete(4,:), ‘LineWidth’, 2)

Title (‘x_4(t)’);

Xlabel(‘t’);

Hl=legend(‘X(t) – реальный сигнал’,’w(t) – возмущающее воздействие’,0);

Set(hl,’FontName’,’Courier’);

Grid on

Figure(10);

Plot(time_X, X(5,:),’r-‘, time_X, w_discrete(5,:), ‘LineWidth’, 2)

Title (‘x_5(t)’);

Xlabel(‘t’);

Hl=legend(‘X(t) – реальный сигнал’,’w(t) – возмущающее воздействие’,0);

Set(hl,’FontName’,’Courier’);

Grid on

Figure(11);

Plot(time_X, q(1,:), time_X, q(2,:), time_X, q(3,:), time_X, q(4,:), time_X, q(5,:), ‘LineWidth’, 2)

Title (‘q(t)- vector-function’);

Xlabel(‘t’);

Hl=legend(‘q_1(t)’, ‘q_2(t)’, ‘q_3(t)’, ‘q_4(t)’, ‘q_5(t)’);

Set(hl,’FontName’,’Courier’);

Grid on

AKOR_slegenie_na_konech_interval_I_podxod. m

Clc

Clear all

Close all

Poryadok = 5;

% ————————————————————————%

B_0 = 5;

B_1 = 9;

% Укороченная система данного объекта

A_5 = 0.1153;

A_4 = 1.78;

A_3 = 3.92;

A_2 = 14.42;

A_1 = 8.583;

A_0 = 0;

% ————————————————————————%

% Приведение системы

B0 = b_0/a_5;

B1 = b_1/a_5;

A5 = a_5/a_5;

A4 = a_4/a_5;

A3 = a_3/a_5;

A2 = a_2/a_5;

A1 = a_1/a_5;

A0 = a_0/a_5;

% ————————————————————————%

% ————————————————————————%

% Представление системы в пространстве состояний

A = [0 1 0 0 0;

0 0 1 0 0;

0 0 0 1 0

0 0 0 0 1;

-a0 – a1 – a2 – a3 – a4];

B = [0; 0; 0; 0; 1];

C = [b0 b1 0 0 0];

% Начальные условия

X_0 = [10; 0; 6; 4; 8;];

Time = 1;

% ————————————————————————%

% ————————————————————————%

% Получение max значений из файла

Load Sostoyaniya X_max U_max

% ————————————————————————%

% Нахождение элементов матриц Q и R

R(1) = 100;

Q(1) = 1/poryadok * r(1) * (U_max)^2 / (X_max(1))^2;

For i = 2 : poryadok

Q(i) = q(1) * (X_max(1))^2 / (X_max(i))^2;

End

Q = diag(q);

R = diag(r);

% Для изменения коэффициентов

% Q(1,1) = Q(1,1)*1e+10;

% Q(2,2) = Q(2,2)*1e+8;

% Q(3,3) = Q(3,3)*1e+6;

% Q(4,4) = Q(4,4)*1e+2;

% Q(5,5) = Q(5,5)*1e+2;

Q(1,1) = Q(1,1)*1e+12;

Q(2,2) = Q(2,2)*1e+8;

Q(3,3) = Q(3,3)*1e+7;

Q(4,4) = Q(4,4)*1e+0;

Q(5,5) = Q(5,5)*1e+2;

R(1,1) = R(1,1);

% ————————————————————————%

% Задающее воздействие

A_o = [0 1 0 0 0;

0 0 1 0 0;

0 0 0 1 0

0 0 0 0 1;

-a0 – a1 – a2 – a3 – a4];

X_o_0 = [12; 10; 14; 8; 16];

% ————————————————————————%

% ————————————————————————%

% Расширенный вектор состояния и расширенные матрицы A, B,Q

%X_rassh = [X_0; X_o];

NULL_M1 = zeros(size(A));

A_rassh = [A NULL_M1;

NULL_M1 A_o];

NULL_M2 = zeros(length(A(:,1)), 1);

B_rassh = [B; NULL_M2];

Q_rassh = [Q – Q;

-Q Q];

X_rassh_0 = [X_0; X_o_0]

% ————————————————————————%

P_nach = zeros(2*poryadok, 2*poryadok);%+ones(poryadok, poryadok);

% ————————————————————————%

% Решение уравнения Риккати методом обратного интегрирования

P = Solve_Riccati_Method_Revers_Integr(A_rassh, B_rassh, Q_rassh, R,Time,2*poryadok, P_nach)

% ————————————————————————%

% Нахождение переменных коэффициентов регулятора

Load Solve_Riccati_Method_Revers_Integr_for_slegenie Time_R P N_str

% ————————————————————————%

% % Формирование матриц P11 и P12

PP = P;

For k = 1 : N_str

P = reshape(PP(k, :), 2*poryadok, 2*poryadok);

For i = 1 : poryadok

For j = 1 : poryadok

P11(i, j,k) = P(i, j);

End

End

For i = 1 : poryadok

For j = (poryadok+1) : (2*poryadok)

P12(i, j-poryadok, k) = P(i, j);

End

End

End

P11(:,:,k)

P12(:,:,k)

% ————————————————————————%

For k = 1 : N_str

K_o(k, 🙂 = – inv(R) * B’ * P11(:,:,k);

K_pr(k, 🙂 = – inv(R) * B’ * P12(:,:,k);

End

% Формирование вектора коэффициентов регулятора

% в прямом порядке

Size(K_o)

Size(K_pr)

I = 1;

Len_K = length(K_o(:,1))

For j = len_K : -1 : 1

K_o_p(i,:) = K_o(j,:)

K_pr_p(i,:) = K_pr(j,:);

I = i + 1;

End

% ————————————————————————%

% Построение графика переменных коэффициентов регулятора обратной связи

% в прямом времени

Figure(2)

Plot(Time_R, K_o(:,1),’-‘,Time_R, K_o(:,2),’-‘,Time_R, K_o(:,3),’-‘,…

Time_R, K_o(:,4),’-‘,Time_R, K_o(:,5),’-‘, ‘LineWidth’, 2);

Xlabel(‘t’)

Tit1 = title(‘Коэффициенты обратной связи в прямом времени’);

Set(tit1,’FontName’,’Courier’);

Hl=legend(‘k_1_о_с’,’k_2_о_с’,’k_3_о_с’,’k_4_о_с’,’k_5_о_с’,0);

Set(hl,’FontName’,’Courier’);

Grid on;

% ————————————————————————%

% ————————————————————————%

% Построение графика переменных коэффициентов регулятора прямой связи

% в прямом времени

Figure(3)

Plot(Time_R, K_pr(:,1),’-‘,Time_R, K_pr(:,2),’-‘,Time_R, K_pr(:,3),’-‘,…

Time_R, K_pr(:,4),’-‘,Time_R, K_pr(:,5),’-‘, ‘LineWidth’, 2);

Xlabel(‘t’)

Tit1 = title(‘Коэффициенты прямой связи в прямом времени’);

Set(tit1,’FontName’,’Courier’);

Hl=legend(‘k_1_п_с’,’k_2_п_с’,’k_3_п_с’,’k_4_п_с’,’k_5_п_с’,0);

Set(hl,’FontName’,’Courier’);

Grid on;

% ————————————————————————%

% ————————————————————————%

% Нахождение отслеживаемого сигнала

X_o(:,1) = X_o_0;

H = 0.01;

For k = 1 : len_K

X_o(:, k+1) = X_o(:, k) + h * A_o * X_o(:, k);

End

X_o(:, k+1) = [];

% ————————————————————————%

% ————————————————————————%

For k = 1 : len_K

A_(:,:,k) = A + B * K_o_p(k,:);

End

Size(A_)

% ————————————————————————%

% Нахождение фазовых координат

X(:,1) = X_0;

Time_X(1) = 0;

For k = 1 : len_K

X(:, k+1) = X(:, k) + h * (A_(:,:,k) * X(:, k) + B * K_pr_p(k,:) * X_o(:,k));

Time_X(k+1) = time_X(k) + h;

End

X(:, k+1) = [];

Time_X(k+1) = [];

% ————————————————————————%

% Нахождение управления

For k = 1 : len_K

U(k) = K_o_p(k,:) * X(:,k) + K_pr_p(k,:) * X_o(:,k);

End

% ————————————————————————%

% Построение u(t) и X(t)

Figure(4);

Plot(time_X, u, ‘r-‘, ‘LineWidth’, 2)

Title (‘u(t)’);

Xlabel(‘t’)

Hl=legend(‘u(t) – управление’,0);

Set(hl,’FontName’,’Courier’);

Grid on

Figure(5);

Plot(time_X, X(1,:),’r-‘, time_X, X_o(1,:), ‘LineWidth’, 2)

Hold on

Title (‘x_1(t)’);

Xlabel(‘t’);

Hl=legend(‘X(t) – слежение’,’X_o(t) – эталон’,0);

Set(hl,’FontName’,’Courier’);

Grid on

Figure(6);

Plot(time_X, X(2,:),’r-‘, time_X, X_o(2,:), ‘LineWidth’, 2)

Title (‘x_2(t)’);

Xlabel(‘t’);

Hl=legend(‘X(t) – слежение’,’X_o(t) – эталон’,0);

Set(hl,’FontName’,’Courier’);

Grid on

Figure(7);

Plot(time_X, X(3,:),’r-‘, time_X, X_o(3,:), ‘LineWidth’, 2)

Title (‘x_3(t)’);

Xlabel(‘t’);

Hl=legend(‘X(t) – слежение’,’X_o(t) – эталон’,0);

Set(hl,’FontName’,’Courier’);

Grid on

Figure(8);

Plot(time_X, X(4,:),’r-‘, time_X, X_o(4,:), ‘LineWidth’, 2)

Title (‘x_4(t)’);

Xlabel(‘t’);

Hl=legend(‘X(t) – слежение’,’X_o(t) – эталон’,0);

Set(hl,’FontName’,’Courier’);

Grid on

Figure(9);

Plot(time_X, X(5,:),’r-‘, time_X, X_o(5,:), ‘LineWidth’, 2)

Title (‘x_5(t)’);

Xlabel(‘t’);

Hl=legend(‘X(t) – слежение’,’X_o(t) – эталон’,0);

Set(hl,’FontName’,’Courier’);

Grid on

AKOR_slegenie_na_konech_interval_II_podxod. m

Clc

Clear all

Close all

Poryadok = 5;

% ————————————————————————%

B_0 = 5;

B_1 = 9;

% Укороченная система данного объекта

A_5 = 0.1153;

A_4 = 1.78;

A_3 = 3.92;

A_2 = 14.42;

A_1 = 8.583;

A_0 = 0;

% ————————————————————————%

% Приведение системы

B0 = b_0/a_5;

B1 = b_1/a_5;

A5 = a_5/a_5;

A4 = a_4/a_5;

A3 = a_3/a_5;

A2 = a_2/a_5;

A1 = a_1/a_5;

A0 = a_0/a_5;

% ————————————————————————%

% ————————————————————————%

% Представление системы в пространстве состояний

A = [0 1 0 0 0;

0 0 1 0 0;

0 0 0 1 0

0 0 0 0 1;

-a0 – a1 – a2 – a3 – a4];

B = [0; 0; 0; 0; 1];

C = [b0 b1 0 0 0];

% Начальные условия

X_0 = [10; 0; 6; 4; 8];

Time = 45;

H = 0.01;

H = 0.8;

% ————————————————————————%

Tic

% ————————————————————————%

% Получение max значений из файла

Load Sostoyaniya X_max U_max

% ————————————————————————%

% Нахождение элементов матриц Q и R

R(1) = 100;

Q(1) = 1/poryadok * r(1) * (U_max)^2 / (X_max(1))^2;

For i = 2 : poryadok

Q(i) = q(1) * (X_max(1))^2 / (X_max(i))^2;

End

Q = diag(q);

R = diag(r);

% Для изменения коэффициентов

% Q(1,1) = Q(1,1)*1e+12;

% Q(2,2) = Q(2,2)*1e+8;

% Q(3,3) = Q(3,3)*1e+7;

% Q(4,4) = Q(4,4)*1e+0;

% Q(5,5) = Q(5,5)*1e+2;

R(1,1) = R(1,1);

% ————————————————————————%

P_nach = zeros(poryadok, poryadok);%+ones(poryadok, poryadok);

% ————————————————————————%

% Решение уравнения Риккати методом обратного интегрирования

P = Solve_Riccati_Method_Revers_Integr(A, B,Q, R,Time, poryadok, P_nach);

Load Solve_Riccati_Method_Revers_Integr_for_slegenie Time_R P N_str

PP = P;

For k = 1 : N_str

P1 = reshape(PP(k, :), poryadok, poryadok);

For i = 1 : poryadok

For j = 1 : poryadok

P2(i, j,k) = P1(i, j);

End

End

End

Size_P = size(P2)

% ————————————————————————%

Tic

% ————————————————————————%

% Получение дискретных значений задающего воздействия в обратном времени

% для нахождения вспомогательной функции q(t)

Zadayushee_Vozdeistvie_Discrete_Revers_Modern(h, 0, Time);

% ————————————————————————%

Load Zadayushee_Vozdeistvie_Discrete_Revers X_o_discrete_rev

% ————————————————————————%

Size(X_o_discrete_rev);

% Нахождение q(t)

For i = 1 : poryadok

Qq = – P_nach(:,:,1) * X_o_discrete_rev(i,1);

Q(i,1) = qq(i,1);

End

% Интегрирование q(t) в обратном времени

For k = 1 : N_str

Q(:, k+1) = q(:, k) – h * ((P2(:,:,k)*B*inv(R)*B’-A’) * q(:, k) + Q*X_o_discrete_rev(:,k));

End

Q(:, k+1) = [];

Size_q = size(q)

% ————————————————————————%

% Нахождение переменных коэффициентов регулятора

For k = 1 : N_str

K_o(k, 🙂 = – inv(R) * B’ * P2(:,:,k);

K_pr(k, 🙂 = – inv(R) * B’;

End

% Формирование вектора коэффициентов регулятора, значений задающего

% воздействия, значений вспомогательной функции в прямом порядке

Size(K_o);

Size(K_pr);

K_pr_p = K_pr;

I = 1;

Len_K = length(K_o(:,1));

For j = len_K : -1 : 1

K_o_p(i,:) = K_o(j,:);

X_o_discrete(:,i) = X_o_discrete_rev(:,j);

Q_pr(:, i) = q(:, j);

I = i + 1;

End

% ————————————————————————%

% Построение графика переменных коэффициентов регулятора обратной связи

% в прямом времени

Toc

Figure(3)

Plot(Time_R, K_o(:,1),’-‘,Time_R, K_o(:,2),’-‘,Time_R, K_o(:,3),’-‘,…

Time_R, K_o(:,4),’-‘,Time_R, K_o(:,5),’-‘, ‘LineWidth’, 2);

Xlabel(‘t’)

Tit1 = title(‘Коэффициенты обратной связи в прямом времени’);

Set(tit1,’FontName’,’Courier’);

Hl=legend(‘k_1_о_с’,’k_2_о_с’,’k_3_о_с’,’k_4_о_с’,’k_5_о_с’,0);

Set(hl,’FontName’,’Courier’);

Grid on;

% ————————————————————————%

% ————————————————————————%

% Построение графика переменных коэффициентов регулятора прямой связи

% в прямом времени

Figure(4)

Plot(Time_R, K_pr(:,1),’-‘,Time_R, K_pr(:,2),’-‘,Time_R, K_pr(:,3),’-‘,…

Time_R, K_pr(:,4),’-‘,Time_R, K_pr(:,5),’-‘, ‘LineWidth’, 2);

Xlabel(‘t’)

Tit1 = title(‘Коэффициенты прямой связи в прямом времени’);

Set(tit1,’FontName’,’Courier’);

Hl=legend(‘k_1_п_с’,’k_2_п_с’,’k_3_п_с’,’k_4_п_с’,’k_5_п_с’,0);

Set(hl,’FontName’,’Courier’);

Grid on;

% ————————————————————————%

Tic

% ————————————————————————%

For k = 1 : len_K

A_(:,:,k) = A + B * K_o_p(k,:);

End

Size_A_ = size(A_)

% ————————————————————————%

% ————————————————————————%

% Нахождение фазовых координат

X(:,1) = X_0;

Time_X(1) = 0;

For k = 1 : len_K

X(:, k+1) = X(:, k) + h * (A_(:,:,k) * X(:, k) + B * K_pr_p(k,:) * q_pr(:,k));

Time_X(k+1) = time_X(k) + h;

End

X(:, k+1) = [];

Time_X(k+1) = [];

Size_X = size(X)

% ————————————————————————%

% ————————————————————————%

% Нахождение управления

For k = 1 : len_K

U(k) = K_o_p(k,:) * X(:,k) + K_pr_p(k,:) * q_pr(:,k);

End

Size_u = size(u)

% ————————————————————————%

Toc

% Построение u(t) и X(t)

Figure(5);

Plot(time_X, u, ‘r-‘, ‘LineWidth’, 2)

Title (‘u(t)’);

Xlabel(‘t’)

Hl=legend(‘u(t) – управление’,0);

Set(hl,’FontName’,’Courier’);

Grid on

Figure(6);

Plot(time_X, X(1,:),’r-‘, time_X, X_o_discrete(1,:), time_X, X_o_discrete(1,:)-0.8,’LineWidth’, 2)

Hold on

Title (‘x_1(t)’);

Xlabel(‘t’);

Hl=legend(‘X(t) – слежение’,’X_o(t) – эталон’, ‘уровень’,0);

Set(hl,’FontName’,’Courier’);

Grid on

Figure(7);

Plot(time_X, X(2,:),’r-‘, time_X, X_o_discrete(2,:), ‘LineWidth’, 2)

Title (‘x_2(t)’);

Xlabel(‘t’);

Hl=legend(‘X(t) – слежение’,’X_o(t) – эталон’,0);

Set(hl,’FontName’,’Courier’);

Grid on

Figure(8);

Plot(time_X, X(3,:),’r-‘, time_X, X_o_discrete(3,:), ‘LineWidth’, 2)

Title (‘x_3(t)’);

Xlabel(‘t’);

Hl=legend(‘X(t) – слежение’,’X_o(t) – эталон’,0);

Set(hl,’FontName’,’Courier’);

Grid on

Figure(9);

Plot(time_X, X(4,:),’r-‘, time_X, X_o_discrete(4,:), ‘LineWidth’, 2)

Title (‘x_4(t)’);

Xlabel(‘t’);

Hl=legend(‘X(t) – слежение’,’X_o(t) – эталон’,0);

Set(hl,’FontName’,’Courier’);

Grid on

Figure(10);

Plot(time_X, X(5,:),’r-‘, time_X, X_o_discrete(5,:), ‘LineWidth’, 2)

Title (‘x_5(t)’);

Xlabel(‘t’);

Hl=legend(‘X(t) – слежение’,’X_o(t) – эталон’,0);

Set(hl,’FontName’,’Courier’);

Grid on

AKOR_slegenie_so_skolz_intervalami_Modern. m

Function AKOR_slegenie_so_skolz_intervalami_Modern

Clc

Clear all

Close all

Poryadok = 5;

% ————————————————————————%

B_0 = 5;

B_1 = 9;

% Укороченная система данного объекта

A_5 = 0.1153;

A_4 = 1.78;

A_3 = 3.92;

A_2 = 14.42;

A_1 = 8.583;

A_0 = 0;

% ————————————————————————%

% Приведение системы

B0 = b_0/a_5;

B1 = b_1/a_5;

A5 = a_5/a_5;

A4 = a_4/a_5;

A3 = a_3/a_5;

A2 = a_2/a_5;

A1 = a_1/a_5;

A0 = a_0/a_5;

% ————————————————————————%

% ————————————————————————%

% Представление системы в пространстве состояний

A = [0 1 0 0 0;

0 0 1 0 0;

0 0 0 1 0

0 0 0 0 1;

-a0 – a1 – a2 – a3 – a4];

B = [0; 0; 0; 0; 1];

C = [b0 b1 0 0 0];

% Начальные условия

X_0 = [10; 0; 6; 4; 8];

Time = 45;

Kolvo_intervalov = 3;

H = 0.01;

H = 0.8;

% ————————————————————————%

% ————————————————————————%

% Получение max значений из файла

Load Sostoyaniya X_max U_max

% ————————————————————————%

% Нахождение элементов матриц Q и R

R(1) = 100;

Q(1) = 1/poryadok * r(1) * (U_max)^2 / (X_max(1))^2;

For i = 2 : poryadok

Q(i) = q(1) * (X_max(1))^2 / (X_max(i))^2;

End

Q = diag(q);

R = diag(r);

% Для изменения коэффициентов

% Q(1,1) = Q(1,1)*1e+13;

% Q(2,2) = Q(2,2)*1e+10;

% Q(3,3) = Q(3,3)*1e+8;

% Q(4,4) = Q(4,4)*1e+5;

% Q(5,5) = Q(5,5)*1e+2;

R(1,1) = R(1,1);

% ————————————————————————%

% ——————Скользящие интервалы———————————-%

Shag = Time/Kolvo_intervalov;

Time1 = shag

Time2 = 2*shag

Time3 = Time

% ————————————————————————%

P_nach = zeros(poryadok, poryadok);%+ones(poryadok, poryadok);

% ————————————————————————%

% Решение уравнения Риккати методом обратного интегрирования

P = Solve_Riccati_Method_Revers_Integr(A, B,Q, R,Time1,poryadok, P_nach);

Load Solve_Riccati_Method_Revers_Integr_for_slegenie Time_R P N_str

PP = P;

For k = 1 : N_str

P1 = reshape(PP(k, :), poryadok, poryadok);

For i = 1 : poryadok

For j = 1 : poryadok

P2(i, j,k) = P1(i, j);

End

End

End

Size_P = size(P2)

% ————————————————————————%

% Нахождение переменных коэффициентов регулятора

For k = 1 : N_str

K_o(k, 🙂 = – inv(R) * B’ * P2(:,:,k);

K_pr(k, 🙂 = – inv(R) * B’;

End

% ————————————————————————%

Tic

% 1 интервал

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Solve_Interval(P_nach, N_str, h, P2, A, B,Q, R, 0, Time1, X_0, poryadok, K_o, K_pr);

Load Solve_Interval time_X X u X_o_discrete

Time_X1 = time_X;

X1 = X;

U1 = u;

X_o_discrete1 = X_o_discrete;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 2 интервал

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Solve_Interval(P_nach, N_str, h, P2, A, B,Q, R, Time1, Time2, X1(:,N_str), poryadok, K_o, K_pr);

Load Solve_Interval time_X X u X_o_discrete

Time_X2 = time_X;

X2 = X;

U2 = u;

X_o_discrete2 = X_o_discrete;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 3 интервал

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Solve_Interval(P_nach, N_str, h, P2, A, B,Q, R, Time2, Time3, X2(:,N_str), poryadok, K_o, K_pr);

Load Solve_Interval time_X X u X_o_discrete

Time_X3 = time_X;

X3 = X;

U3 = u;

X_o_discrete3 = X_o_discrete;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Toc

% ————————————————————————%

% Объединение интервалов

Time_X = [time_X1 time_X2 time_X3];

U = [u1 u2 u3];

X = [X1 X2 X3];

X_o_discrete = [X_o_discrete1 X_o_discrete2 X_o_discrete3];

% ————————————————————————%

% ————————————————————————%

% Построение u(t) и X(t)

Figure(3);

Plot(time_X, u, ‘r-‘,’LineWidth’, 2);

Title (‘u(t)’);

Xlabel(‘t’)

Hl=legend(‘u(t) – управление’,0);

Set(hl,’FontName’,’Courier’);

Grid on

Figure(4);

Plot(time_X, X(1,:),’r-‘, time_X, X_o_discrete(1,:), time_X, X_o_discrete(1,:)-0.8,’LineWidth’, 2)

Hold on

Title (‘x_1(t)’);

Xlabel(‘t’);

Hl=legend(‘X(t) – слежение’,’X_o(t) – эталон’, ‘уровень’,0);

Set(hl,’FontName’,’Courier’);

Grid on

Figure(5);

Plot(time_X, X(2,:),’r-‘, time_X, X_o_discrete(2,:), ‘LineWidth’, 2)

Title (‘x_2(t)’);

Xlabel(‘t’);

Hl=legend(‘X(t) – слежение’,’X_o(t) – эталон’,0);

Set(hl,’FontName’,’Courier’);

Grid on

Figure(6);

Plot(time_X, X(3,:),’r-‘, time_X, X_o_discrete(3,:), ‘LineWidth’, 2)

Title (‘x_3(t)’);

Xlabel(‘t’);

Hl=legend(‘X(t) – слежение’,’X_o(t) – эталон’,0);

Set(hl,’FontName’,’Courier’);

Grid on

Figure(7);

Plot(time_X, X(4,:),’r-‘, time_X, X_o_discrete(4,:), ‘LineWidth’, 2)

Title (‘x_4(t)’);

Xlabel(‘t’);

Hl=legend(‘X(t) – слежение’,’X_o(t) – эталон’,0);

Set(hl,’FontName’,’Courier’);

Grid on

Figure(8);

Plot(time_X, X(5,:),’r-‘, time_X, X_o_discrete(5,:), ‘LineWidth’, 2)

Title (‘x_5(t)’);

Xlabel(‘t’);

Hl=legend(‘X(t) – слежение’,’X_o(t) – эталон’,0);

Set(hl,’FontName’,’Courier’);

Grid on

Function Solve_Interval(P_nach, N_str, h, P2, A, B,Q, R, T_nach, T_konech, X_0, poryadok, K_o, K_pr)

Zadayushee_Vozdeistvie_Discrete_Revers_Modern(h, T_nach, T_konech);

Load Zadayushee_Vozdeistvie_Discrete_Revers X_o_discrete_rev

% ————————————————————————%

% Нахождение q(t)

For i = 1 : poryadok

Qq = – P_nach(:,:,1) * X_o_discrete_rev(i,1);

Q(i,1) = qq(i,1);

End

% Интегрирование q(t) в обратном времени

For k = 1 : N_str

Q(:, k+1) = q(:, k) – h * ((P2(:,:,k)*B*inv(R)*B’-A’) * q(:, k) + Q*X_o_discrete_rev(:,k));

End

Q(:, k+1) = [];

Size_q = size(q)

% ————————————————————————%

% Формирование вектора коэффициентов регулятора, значений задающего

% воздействия, значений вспомогательной функции в прямом порядке

K_pr_p = K_pr;

I = 1;

For j = N_str : -1 : 1

K_o_p(i,:) = K_o(j,:);

X_o_discrete(:,i) = X_o_discrete_rev(:,j);

Q_pr(:, i) = q(:, j);

I = i + 1;

End

% ————————————————————————%

% ————————————————————————%

For k = 1 : N_str

A_(:,:,k) = A + B * K_o_p(k,:);

End

Size_A_ = size(A_)

% ————————————————————————%

% ————————————————————————%

% Нахождение фазовых координат

X(:,1) = X_0;

Time_X(1) = T_nach;

For k = 1 : N_str

X(:, k+1) = X(:, k) + h * (A_(:,:,k) * X(:, k) + B * K_pr_p(k,:) * q_pr(:,k));

Time_X(k+1) = time_X(k) + h;

End

X(:, k+1) = [];

Time_X(k+1) = [];

Size_X = size(X)

% ————————————————————————%

% ————————————————————————%

% Нахождение управления

For k = 1 : N_str

U(k) = K_o_p(k,:) * X(:,k) + K_pr_p(k,:) * q_pr(:,k);

End

Size_u = size(u)

Save Solve_Interval time_X X u X_o_discrete

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Sintez_nablyud_polnogo_poryadka. m

Clc

Clear all

Close all

Poryadok = 5;

% ————————————————————————%

B_0 = 5;

B_1 = 9;

% Укороченная система данного объекта

A_5 = 0.1153;

A_4 = 1.78;

A_3 = 3.92;

A_2 = 14.42;

A_1 = 8.583;

A_0 = 0;

% ————————————————————————%

% Приведение системы

B0 = b_0/a_5;

B1 = b_1/a_5;

A5 = a_5/a_5;

A4 = a_4/a_5;

A3 = a_3/a_5;

A2 = a_2/a_5;

A1 = a_1/a_5;

A0 = a_0/a_5;

% ————————————————————————%

% ————————————————————————%

% Представление системы в пространстве состояний

A = [0 1 0 0 0;

0 0 1 0 0;

0 0 0 1 0;

0 0 0 0 1;

-a0 – a1 – a2 – a3 – a4]

B = [0; 0; 0; 0; 1]

C = [b0 b1 0 0 0]

% Начальные условия

X_0 = [10; 0; 6; 4; 8]

Time = 10;

% ————————————————————————%

% ————————————————————————%

% Получение max значений из файла

Load Sostoyaniya X_max U_max

% ————————————————————————%

% Нахождение элементов матриц Q и R

R(1) = 100;

Q(1) = 1/poryadok * r(1) * (U_max)^2 / (X_max(1))^2;

For i = 2 : poryadok

Q(i) = q(1) * (X_max(1))^2 / (X_max(i))^2;

End

Q = diag(q)

R = diag(r)

% Для изменения коэффициентов

Q(1,1) = Q(1,1);

Q(2,2) = Q(2,2);

Q(3,3) = Q(3,3);

Q(4,4) = Q(4,4);

Q(5,5) = Q(5,5);

% Q(1,1) = Q(1,1)*1e+12;

% Q(2,2) = Q(2,2)*1e+8;

% Q(3,3) = Q(3,3)*1e+7;

% Q(4,4) = Q(4,4)*1e+0;

% Q(5,5) = Q(5,5)*1e+2;

R(1,1) = R(1,1);

% ————————————————————————%

P_nach = zeros(poryadok, poryadok);%+ones(poryadok, poryadok);

% ————————————————————————%

% Решение уравнения Риккати методом обратного интегрирования

P1 = Solve_Riccati_Method_Revers_Integr(A, B,Q, R,Time, poryadok, P_nach)

% ————————————————————————%

% Построение графика коэффициентов регулятора

Load Solve_Riccati_Method_Revers_Integr Time_R P N_str

PP = P;

For i = 1 : N_str

P = reshape(PP(i, :), poryadok, poryadok);

K(i, 🙂 = – inv(R)*B’*P;

End

Figure(2)

Plot(Time_R, K(:,1),’-‘,Time_R, K(:,2),’-‘,Time_R, K(:,3),’-‘,Time_R, K(:,4),’-‘,Time_R, K(:,5),’-‘, ‘LineWidth’, 2);

Xlabel(‘t’)

Tit1 = title(‘Коэффициенты обратной связи в прямом времени’);

Set(tit1,’FontName’,’Courier’);

Hl=legend(‘k_1_о_с’,’k_2_о_с’,’k_3_о_с’,’k_4_о_с’,’k_5_о_с’,0);

Set(hl,’FontName’,’Courier’);

Grid on;

% ————————————————————————%

% Нахождение коэффициентов регулятора

Disp(‘Коэффициенты регулятора:’)

K = – inv(R) * B’ * P1

% ————————————————————————%

% ————————————————————————%

A_ = A + B * K;

% ————————————————————————%

% Нахождение фазовых координат

X(:,1) = X_0;

H = 0.01;

Time_X(1) = 0;

For k = 1 : N_str

X(:, k+1) = X(:, k) + h * A_ * X(:, k);

Time_X(k+1) = time_X(k) + h;

End

X(:, k+1) = [];

Time_X(k+1) = [];

% ————————————————————————%

% Нахождение управления

For k = 1 : N_str

U(k) = K * X(:,k);

End

% ————————————————————————%

% Нахождение коэффициентов наблюдателя

M_n = [C’ A’*C’ (A^2)’*C’ (A^3)’*C’ (A^4)’*C’]

Rank_M_n = rank(M_n)

A_r = A_

Disp(‘Спектр матрицы регулятора:’)

Spektr_A_r = eig(A_r)

Koeff = 1;

Min_lyamda_A_r = min(real(spektr_A_r))

% lyamda = min_lyamda_A_r * koeff;

Lyamda = -5;

Disp(‘Спектр матрицы наблюдателя эталонный:’)

Lyamda_A_n = [lyamda – koeff * 4; lyamda – koeff * 3; lyamda – koeff * 2;…

Lyamda – koeff; lyamda]’

Syms k_n1 k_n2 k_n3 k_n4 k_n5 lyam

K_n = [k_n1; k_n2; k_n3; k_n4; k_n5];

Koeff_poly_n_etalon = poly(lyamda_A_n)

Disp(‘Характеристический полином наблюдателя эталонный:’)

Poly_n_etalon = poly2sym(Koeff_poly_n_etalon, lyam)

Disp(‘Характеристический полином наблюдателя реальный:’)

Poly_n_real = collect(expand(simplify(det(lyam*eye(poryadok) – (A – K_n*C)))),lyam)

Raznost_poly = collect(poly_n_etalon-poly_n_real, lyam)

For i = 1 : poryadok

Koeff_raznost_poly(i) = subs(diff(raznost_poly, poryadok-i, lyam)/factorial(poryadok-i),lyam,0);

End

Koeff_raznost_poly

[Kn1 Kn2 Kn3 Kn4 Kn5]= solve(Koeff_raznost_poly(5), Koeff_raznost_poly(4),…

Koeff_raznost_poly(3), Koeff_raznost_poly(2), Koeff_raznost_poly(1), …

K_n1, k_n2, k_n3, k_n4, k_n5)

Kn = [Kn1; Kn2; Kn3; Kn4; Kn5];

Kn = vpa(Kn,50)

% Проверка

Proverka = solve(det(lyam*eye(poryadok)-(A-Kn*C)))

% ————————————————————————%

% ————————————————————————%

% Нахождение x и x_оценочного

X_ocen_0 = [0 0 0 0 0]’;

A_rash = [A B*K;

Kn*C A-Kn*C+B*K]

X_rash_0 = [X_0;X_ocen_0]

X_rash(:,1) = X_rash_0;

For k = 1 : N_str

X_rash(:,k+1) = X_rash(:,k) + h * A_rash * X_rash(:,k);

End

X_rash(:,k+1) = [];

% ————————————————————————%

% ————————————————————————%

% Разделение x и x_оценочного

For i = 1 : poryadok

X_n(i,:) = X_rash(i,:);

End

For i = poryadok + 1 : 2*poryadok

X_n_ocen(i – poryadok,:) = X_rash(i,:);

End

% ————————————————————————%

% ————————————————————————%

% Нахождение управления

For i = 1 : N_str

U_n(i) = K * X_n_ocen(:,i);

End

% Построение u(t) и X(t)

Figure(3);

Plot(time_X, u, ‘r-‘, time_X, u_n, ‘b-‘, ‘LineWidth’, 2)

Title (‘u(t)’);

Xlabel(‘t’)

Hl=legend(‘управление без наблюдателя’,’управление c наблюдателем’);

Set(hl,’FontName’,’Courier’);

Grid on

Figure(4);

Plot(time_X, X(1,:), time_X, X_n(1,:), time_X, X_n_ocen(1,:),’LineWidth’, 2)

Hold on

Title (‘x_1(t)’);

Xlabel(‘t’)

Hl=legend(‘x_1(t) без наблюдателя’,’x_1(t) c наблюдателем’, ‘x_о_ц_е_н_1(t)’);

Set(hl,’FontName’,’Courier’);

Grid on

Figure(5);

Plot(time_X, X(2,:), time_X, X_n(2,:), time_X, X_n_ocen(2,:),’LineWidth’, 2)

Title (‘x_2(t)’);

Xlabel(‘t’)

Hl=legend(‘x_2(t) без наблюдателя’,’x_2(t) c наблюдателем’, ‘x_о_ц_е_н_2(t)’);

Set(hl,’FontName’,’Courier’);

Grid on

Figure(6);

Plot(time_X, X(3,:), time_X, X_n(3,:), time_X, X_n_ocen(3,:),’LineWidth’, 2)

Title (‘x_3(t)’);

Xlabel(‘t’)

Hl=legend(‘x_3(t) без наблюдателя’,’x_3(t) c наблюдателем’, ‘x_о_ц_е_н_3(t)’);

Set(hl,’FontName’,’Courier’);

Grid on

Figure(7);

Plot(time_X, X(4,:), time_X, X_n(4,:), time_X, X_n_ocen(4,:),’LineWidth’, 2)

Title (‘x_4(t)’);

Xlabel(‘t’)

Hl=legend(‘x_4(t) без наблюдателя’,’x_4(t) c наблюдателем’, ‘x_о_ц_е_н_4(t)’);

Set(hl,’FontName’,’Courier’);

Grid on

Figure(8);

Plot(time_X, X(5,:), time_X, X_n(5,:), time_X, X_n_ocen(5,:),’LineWidth’, 2)

Title (‘x_5(t)’);

Xlabel(‘t’)

Hl=legend(‘x_5(t) без наблюдателя’,’x_5(t) c наблюдателем’, ‘x_о_ц_е_н_5(t)’);

Set(hl,’FontName’,’Courier’);

Grid on

Solve_Riccati_Method_Diag. m

% ————————————————————————%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Метод диагонализации для решения алгебраического уравнения Риккати

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Function P = Solve_Riccati_Method_Diag(A, B,Q, R)

% Расширенная матрица системы

Z = [A B*inv(R)*B’;

Q – A’]

% Нахождение собственных векторов и собственных чисел матрицы Z

[V, D] = eig(Z)

% ————————————————————————%

% Построение матрицы S

% Индексы столбцов собственных значений Re(lyamda) > 0

Ind_Re_plus = find(sum(real(D)) > 0);

% Индексы столбцов собственных значений Re(lyamda) < 0

Ind_Re_minus = find(sum(real(D)) < 0);

% Формирование матрицы D в виде Re(lyamda) > 0 -> Re(lyamda) < 0

D1 = sum(D(:, Ind_Re_plus));

D2 = sum(D(:, Ind_Re_minus));

D = [D1 D2];

% Формирование матрицы S в виде Re(lyamda) > 0 -> Re(lyamda) < 0

S1 = V(:, Ind_Re_plus);

S2 = V(:, Ind_Re_minus);

S = [S1 S2];

% Поиск столбцов с комплексными корнями в матрице D

Ind_Complex_D = find(imag(D) ~= 0);

% Формирование конечной матрицы S

For i = 1 : 2 : length(Ind_Complex_D)

S (:, Ind_Complex_D(i) + 1) = imag(S(:, Ind_Complex_D(i)));

S (:, Ind_Complex_D(i)) = real(S(:, Ind_Complex_D(i)));

End

S = S

% ————————————————————————%

Poryadok = length(A(1,:));

S12 = S(1 : poryadok, poryadok+1 : 2*poryadok);

S22 = S(poryadok+1 : 2*poryadok, poryadok+1 : 2*poryadok);

% ————————————————————————%

% Вычисление матрицы P

P = – S22 * inv(S12);

Solve_Riccati_Method_Revers_Integr. m

% ————————————————————————%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Решение уравнения Риккати интегрированием в обратном времени

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Function P = Solve_Riccati_Method_Revers_Integr(A, B,Q, R,Time, poryadok, P1)

Save For_Riccati A B Q R poryadok

% Решение дифференциального уравнения Риккати

P1 = reshape(P1, poryadok^2, 1);

[Time_R, P] = ode45(@Riccati, [Time : -0.01 : 0], P1);

[N_str, N_stolb] = size(P);

% Построение полученного решения

Figure(1)

For i = 1 : poryadok^2

Plot(Time_R, P(:,i),’-‘)

Hold on

End

% plot(Time_R, P(:,1),’-‘,Time_R, P(:,2),’-‘,Time_R, P(:,3),’-‘,Time_R, P(:,4),’-‘,Time_R, P(:,5),’-‘,Time_R, P(:,6),’-‘,…

% Time_R, P(:,7),’-‘,Time_R, P(:,8),’-‘,Time_R, P(:,9),’-‘,Time_R, P(:,10),’-‘,Time_R, P(:,11),’-‘,Time_R, P(:,12),’-‘,…

% Time_R, P(:,13),’-‘,Time_R, P(:,14),’-‘,Time_R, P(:,15),’-‘,Time_R, P(:,16),’-‘,Time_R, P(:,17),’-‘,Time_R, P(:,18),’-‘,…

% Time_R, P(:,19),’-‘,Time_R, P(:,20),’-‘,Time_R, P(:,21),’-‘,Time_R, P(:,22),’-‘,Time_R, P(:,23),’-‘,Time_R, P(:,24),’-‘,…

% Time_R, P(:,25),’-‘, ‘lineWidth’, 2);

Grid on;

Tit1 = title(‘Решения уравнения Риккати’);

Set(tit1,’FontName’,’Courier’);

Xlabel(‘t’);

% legend(‘p_1′,’p_2′,’p_3′,’p_4′,’p_5′,’p_6′,’p_7′,’p_8′,’p_9′,’p_1_0′,’p_1_1′,’p_1_2′,’p_1_3′,’p_1_4′,’p_1_5′,’p_1_6’,…

% ‘p_1_7′,’p_1_8′,’p_1_9′,’p_2_0′,’p_2_1′,’p_2_2′,’p_2_3′,’p_2_4′,’p_2_5’);

Save Solve_Riccati_Method_Revers_Integr Time_R P N_str

Save Solve_Riccati_Method_Revers_Integr_for_slegenie Time_R P N_str

P = reshape(P(N_str,:), poryadok, poryadok);

Function dP = Riccati(Time, P)

Load For_Riccati A B Q R poryadok

P = reshape(P, poryadok, poryadok);

% Дифференциальное уравнение Риккати

DP = – P*A – A’*P + P*B*inv(R)*B’*P – Q;

DP = reshape(dP, poryadok^2, 1);

Vozmyshyayushee_Vozdeistvie_Discrete_Revers. m

% Получение дискретных значений возмущающего воздействия в обратном времени

% для нахождения вспомогательной функции q(t)

Function Vozmyshyayushee_Vozdeistvie_Discrete_Revers(h, T_nach, T_konech)

% ————————————————————————%

% Возмущающее воздействие

A = 1;

W = 4*pi;

K = 1;

RETURN = 1;

While RETURN == 1

Disp(‘Возмущающее воздействие – const: 1’)

Disp(‘Возмущающее воздействие – A*sin(w*t): 2’)

Reply = input(‘Выберете возмущающее воздействие [1 или 2]: ‘, ‘s’);

Switch reply

case ‘1’

disp(‘Возмущающее воздействие – const’)

for t = T_konech: – h : T_nach

W_discrete_rev(:, k) = [A + 0 * t; 0; 0; 0; 0];

K = k + 1;

end

RETURN = 2;

case ‘2’

disp(‘Возмущающее воздействие – A*sin(w*t)’)

for t = T_konech: – h : T_nach

W_discrete_rev(:, k) = [A * sin(w * t); 0; 0; 0; 0];

K = k + 1;

end

RETURN = 2;

otherwise

disp(‘Неизвестное воздействие.’)

RETURN = 1;

End

End

Figure(2)

T = T_konech : – h : T_nach;

Plot(t, w_discrete_rev(1,:), ‘r-‘, ‘LineWidth’, 2);

Xlabel(‘t’)

Tit1 = title(‘Возмущающее воздействие’);

Set(tit1,’FontName’,’Courier’);

Hl=legend(‘Возмущающее воздействие’,0);

Set(hl,’FontName’,’Courier’);

Grid on;

Save Vozmyshyayushee_Vozdeistvie_Discrete_Revers w_discrete_rev

% ————————————————————————%

Zadayushee_Vozdeistvie_Discrete_Revers_Modern. m

% Получение дискретных значений задающего воздействия в обратном времени

% для нахождения вспомогательной функции q(t)

Function Zadayushee_Vozdeistvie_Discrete_Revers_Modern(h, T_nach, T_konech)

% ————————————————————————%

% Задающее воздействие

Alfa = 0.2;

Beta = 10;

H = 0.8;

K = 1;

For t = T_konech: – h : T_nach

X_o_1 = 10*exp(-1/5*t)*t+4/5;

X_o_2 = -2*exp(-1/5*t)*t+10*exp(-1/5*t);

X_o_3 = 2/5*exp(-1/5*t)*t-4*exp(-1/5*t);

X_o_4 = -2/25*exp(-1/5*t)*t+6/5*exp(-1/5*t);

X_o_5 = 2/125*exp(-1/5*t)*t-8/25*exp(-1/5*t);

X_o_discrete_rev(:, k) = [X_o_1; X_o_2; X_o_3; X_o_4; X_o_5];

k = k + 1;

End

Figure(2)

T = T_konech : – h : T_nach;

Plot(t, X_o_discrete_rev(1,:), ‘r-‘, t, X_o_discrete_rev(1,:)-H, ‘LineWidth’, 2);

Xlabel(‘t’)

Tit1 = title(‘Задающее воздействие’);

Set(tit1,’FontName’,’Courier’);

Hl=legend(‘Отслеживание зад. возд. на H ‘,’Задающее воздействие’,0);

Set(hl,’FontName’,’Courier’);

Grid on;

Save Zadayushee_Vozdeistvie_Discrete_Revers X_o_discrete_rev

% ————————————————————————%


Зараз ви читаєте: Математическая модель в пространстве состояний линейного стационарного объекта управления