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

КурсоваяПомощь в написанииУзнать стоимостьмоей работы

U2 = -20.605 579 750 692 850 612 199 384 678 400*exp (-40.749 492 463 732 569 497 044 913 225 728+13.583 164 154 577 523 148 853 597 962 240*t)+19.11 167 813 350 479 568 880 577 544 192*exp (-2.544 534 472 800 777 276 275 852 050 432+.68 481 781 576 002 594 545 273 472 548 864*t)+1.3 356 706 538 317 879 411 306 290 741 248*exp… Читать ещё >

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

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 Матрица Фробениуса

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

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

(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-го знака после запятой. Результаты, полученные при переводе системы из начальной точки в конечную, полностью совпадают.

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

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

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

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

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

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

или

Ниже представлен график оптимального управления полученного с помощью скрипта 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⁴ + a3 * s³ + a2 * s² + 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² + b1².*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² + b1².*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² + b1².*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²)*(s-s5))), s)

%pretty (W_s_etal)

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

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

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

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

Chislit_W_s =collect (Slag₁ + Slag₂ + Slag₃ + Slag₄, s);

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

MS =double ([1 1 1 1 0;

6 753 029 497/515578134 -513 659/1058682 10 560 977/850789 4 210 795/295122 1;

77 456 808 434 995 514 797 195 264/126764366837761533378822144 1 874 500 571 398 143 916 113 920/260296441145300889894912 -3 300 780 600 401 724 994 224 128/418364246989311991349248 915 075/98374 4 210 795/295122;

26 189 071 674 868 424 525 008 076 800/253528733675523066757644288 2 853 037 197 681 682 523 619 328/520592882290601779789824 45 476 725 452 203 198 204 870 656/418364246989311991349248 0 915 075/98374;

6 290 947 020 888 110 049 943 093 248/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₁ = simplify (collect (Koeff (1)*(s-s2)*((s+alfa)^2+beta²)*(s-s5), s));

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

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

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

Chislit_W_s =collect ((Slag₁ + Slag₂ + Slag₃ + Slag₄), s);

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

W_s = collect (simplify (Koeff (1)/(s-s1)+Koeff (2)/(s-s2)+(Koeff (4)*s+Koeff (5))/((s+alfa)^2+beta²)+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²+beta²) -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₀ = [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₀;

% Вычисление произведения 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

% Построение 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

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

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

% Формирование 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

[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

%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

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

X = zeros (n, N);

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

for i = 2: N

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

end

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

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

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

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

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

XX = [X₀ X];

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

normaXX = norm (XX (, N))

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

R = abs (X_N — F^N * X₀) — FG * u';

R = R';

z = sum®;

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

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

t_perevoda = N * h;

end

N_opt

t_perevoda

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

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

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

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

figure (1)

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

plot (t, X1, 'b', 'LineWidth', 2);

hl=legend ('x₁(t)');

set (hl, 'FontName', 'Courier');

xlabel ('t, cek'); ylabel ('x₁(t)');

grid on

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

figure (2)

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

plot (t, X2, 'b', 'LineWidth', 2);

hl=legend ('x₂(t)');

set (hl, 'FontName', 'Courier');

xlabel ('t, cek'); ylabel ('x₂(t)');

grid on

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

figure (3)

plot (X1, X2, 'm', 'LineWidth', 2);

hl=legend ('x₂ = x₂(x₁)');

set (hl, 'FontName', 'Courier');

xlabel ('x₁(t)'); ylabel ('x₂(x₁(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₀ = 5;

b₁ = 9;

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

a₅ = 0.1153;

a₄ = 1.78;

a₃ = 3.92;

a₂ = 14.42;

a₁ = 8.583;

a₀ = 0;

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

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

b0 = b₀/a₅;

b1 = b₁/a₅;

a5 = a₅/a₅;

a4 = a₄/a₅;

a3 = a₃/a₅;

a2 = a₂/a₅;

a1 = a₁/a₅;

a0 = a₀/a₅;

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

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

poryadok = 5;

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

Y₀ = [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₀ = B_' * inv (B_ * B_') * Y₀

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⁵ + a4*s⁴ + a3*s³ + a2*s² + a1*s + a0);

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

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

% ИПФ

K_t = simplify (vpa (ilaplace (1 / (a5*s⁵ + a4*s⁴ + a3*s³ + …

a2*s² + 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₀ = double (subs (Matrix_Vron, t,0));

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

T = 3;

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

% vpa (Matrix_Vron_t₀, 6)

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

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

C_ = inv (Matrix_Vron_t₀) * X₀;

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

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

K_Tt₁ = subs (K_t, t, T — t);

K_Tt = diff (K_t);

K_Tt₂ = subs (K_Tt, t, T — t);

K_Ttt = diff (K_Tt);

K_Tt₃ = subs (K_Ttt, t, T — t);

K_Tttt = diff (K_Ttt);

K_Tt₄ = subs (K_Tttt, t, T — t);

K_Ttttt = diff (K_Tttt);

K_Tt₅ = subs (K_Ttttt, t, T — t);

h1_Tt = K_Tt₁

h2_Tt = K_Tt₂

h3_Tt = K_Tt₃

h4_Tt = K_Tt₄

h5_Tt = K_Tt₅

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

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

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₀;

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

% 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₂ = 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₂ = vpa (expand (subs (d_v₂, ks1)), 50);

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

eq₁= diff (d_v₂, ks2);

eq₂= diff (d_v₂, ks3);

eq₃= diff (d_v₂, ks4);

eq₄= diff (d_v₂, ks5);

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

ksi= solve (eq₁, eq₂, eq₃, eq₄);

% Полученные значения 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₂ = 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₂ = double (d_v₂)

gamma_v₂ = 1/d_v₂

% gamma_v₂ = double (gamma_v₂)

u = vpa (expand (simplify (gamma_v₂ * (ks1*h1_Tt + ks2*h2_Tt + ks3*h3_Tt + …

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

% u = vpa (u, 6)

u₀ = 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₀;

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₀ = 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₀ = double (vpa (X₀ — X_real₀, 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₀ = 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₀ = double (vpa (Y₀ — Y_real₀, 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₁(t)');

grid on

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

title ('x₂(t)');

grid on

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

title ('x₃(t)');

grid on

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

title ('x₄(t)');

grid on

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

title ('x₅(t)');

grid on

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

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

title ('y₁(t)');

grid on

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

title ('y₂(t)');

grid on

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

title ('y₃(t)');

grid on

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

title ('y₄(t)');

grid on

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

Gramian_Uprav.m

clc

close all

clear all

format long

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

b₀ = 5;

b₁ = 9;

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

a₅ = 0.1153;

a₄ = 1.78;

a₃ = 3.92;

a₂ = 14.42;

a₁ = 8.583;

a₀ = 0;

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

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

b0 = b₀/a₅;

b1 = b₁/a₅;

a5 = a₅/a₅;

a4 = a₄/a₅;

a3 = a₃/a₅;

a2 = a₂/a₅;

a1 = a₁/a₅;

a0 = a₀/a₅;

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

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

poryadok = 5;

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

Y₀ = [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₀ = B_' * inv (B_ * B_') * Y₀

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²*B A³*B A⁴*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₀))), 50)

u₀ = 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.605 579 750 692 850 612 199 384 678 400*exp (-40.749 492 463 732 569 497 044 913 225 728+13.583 164 154 577 523 148 853 597 962 240*t)+19.11 167 813 350 479 568 880 577 544 192*exp (-2.544 534 472 800 777 276 275 852 050 432+.68 481 781 576 002 594 545 273 472 548 864*t)+1.3 356 706 538 317 879 411 306 290 741 248*exp (-1.7 550 088 311 372 149 728 872 572 649 472+.58 500 294 371 240 501 635 562 581 524 480*t)*cos (-8.3 032 397 968 812 277 195 784 104 968 192+2.7 677 465 989 604 092 531 158 867 247 104*t)+7.2 830 359 327 562 658 351 629 692 043 264*exp (-1.7 550 088 311 372 149 728 872 572 649 472+.58 500 294 371 240 501 635 562 581 524 480*t)*sin (-8.3 032 397 968 812 277 195 784 104 968 192+2.7 677 465 989 604 092 531 158 867 247 104*t)-8.6 096 491 449 877 801 531 228 326 723 584;

Показать весь текст

Заполнить форму текущей работой