Чисельні методи інтегрування та оптимізації складних систем

Московський державний технічний університет ім. Н.Е. Баумана

Калузький філія

Кафедра "САУ та Електротехніки"

ЕІУ3-КФ

Розрахунково-пояснювальна записка до курсової роботи

на тему:

"Чисельні методи інтегрування та оптимізації складних систем"

за курсом:

Системи аналітичних обчислень

Калуга 2007

Федеральне агентство з освіти

Державна освітня установа вищої професійної освіти

«Московський державний технічний університет

імені Н.Е. Баумана »

Калузький філія

Факультет електроніки, інформатики та управління

Кафедра "Системи автоматичного управління та електротехніка" (ЕІУ3-КФ)

З А Д А Н Н Я

на курсовий проект (роботу)

за курсом _____Сістеми аналітичних вичісленій____________

Студент ________Герасімов Е.І._______ група ___САУ-62_________

(Прізвище, ініціали)

Руководитель_________________Корнюшин Ю.П.____________________

(Прізвище, ініціали)

Тема проекту (роботи) Чисельні методи інтегрування і______ оптимізації складних систем

Технічне завдання

Завдання 1.

Практичне вивчення чисельних методів розв'язання нелінійних рівнянь (метод простих ітерацій) і рішення заданого рівняння третього порядку з метою дослідження стійкості заданої системи.

Завдання 2.

Побудова годографа АФЧХ, графіків АЧХ і ФЧХ із зазначенням частот.

Завдання 3.

Практичне вивчення чисельних методів інтегрування диференціальних рівнянь високого порядку (метод Рунге-Кутта 5-го порядку, неявний метод Адамса 4-го порядку) і побудова перехідних процесів.

Завдання 4.

Проведення аналізу заданої системи з використанням спектрального методу (базис: поліноми Чебишева 2-го роду).

Завдання 5.

Практичне вивчення чисельних методів оптимізації (метод Хука-Дживса з використанням методу Фібоначчі) та визначення параметрів коригувального устрою, шляхом мінімізації функціонала якості.

Обсяг і зміст проекту (роботи)

Графічні роботи на ___5_____ аркушах формату ___A3____

Розрахунково-пояснювальна записка на __53____ аркушах формату А4

Структура розрахунково-пояснювальної записки

Обкладинка, Завдання, Зміст, Вступ, Основна частина, Висновок, Література, Додаток (я).

Зміст і структура Основний частини визначається студентом за погодженням з керівником.

Малюнки, таблиці, література оформляються у відповідності з ГОСТ 2.105-89 ЕСКД. Загальні вимоги до текстових документів, ГОСТ 7.32-90 Звіт про науково-дослідній роботі. Загальні вимоги та правила оформлення.

Рекомендована література

Н.Д. Єгупов, Ю.П. Корнюшин, Ю.Л. Лукашенко, А.А. Самохвалов, М.М. Чайковський Складні системи автоматичного управління зі змінними параметрами: алгоритмічне і програмне забезпечення вирішення завдань дослідження та синтезу, Калуга, 2003

Вержбицький. Чисельні методи.

Методи класичної та сучасної теорії автоматичного керування: Підручник у 5-ти т., 2-е вид., Перераб. і доп. Т.3: Синтез регуляторів систем автоматичного управління / За редакцією К.А. Пупкова та Н.Д. Єгупова. - М.: Видавництво МГТУ ім. Н.Е. Баумана, 2004. - 616с.; Іл.

Конспект лекцій з курсу "Системи аналітичних обчислень" за I і II семестр.

Керівник проекту ____________________________

підпис

"______" ________________ 2007

Студент ____________________________

Підпис

"____" _________________ 2007

Зміст

1. Постановка завдання

АНАЛІЗ

Чисельні методи інтегрування

(Дослідження стійкості САУ)

Для заданої системи потрібно визначити:

Передавальну функцію замкненої системи, для випадку ;

Коріння характеристичного рівняння, використовуючи метод січних;

Знайти аналітичні вирази для АЧХ, ФЧХ, АФЧХ;

Побудувати годограф АФЧХ і графіки АЧХ і ФЧХ із зазначенням частот;

Отримати ДУ, що описує дану систему;

Уявити ДУ в нормальній формі Коші;

Знайти аналітичне рішення ДУ;

Знайти рішення ДУ чисельним методом (метод Рунге-Кутта 5-го порядку і метод Адамса неявний 4-го порядку);

Аналіз заданої системи з використанням спектрального методу (базис: поліноми Чебишева 2-го роду).

СИНТЕЗ

Чисельні методи оптимізації

Записати передавальну функцію замкненої системи, з урахуванням того що ;

Отримати ДУ, що описує дану систему;

Уявити ДУ в нормальній формі Коші;

Обчислити невідомі параметри коригувального пристрою мінімізуючи функціонал якості виду методом Хука-Дживса з використанням методу Фібоначчі. Для знаходження реальної передавальної характеристики системи необхідно використовувати один з методів чисельного інтегрування.

Провести аналіз отриманих результатів.

Визначити невідомі параметри коригувального пристрою , Що забезпечує Робастне якість сімейства систем.

АНАЛІЗ

Вихідні дані: структурна схема заданої системи зображена на рис. 1, а значення параметрів системи наведені в таблиці 1.

у (t) x (t)

- -

Рис. 1. Структурна схема системи.

Таблиця 1

K1	K2	K3	T1, c	T2, c	T3, c
15	10	1	1.2	0.3	0.7

2.1 Передавальна функція замкнутої системи, для випадку

Передавальною функцією (ПФ) САУ називається відношення перетворення Лапласа сигналу на виході системи до перетворення Лапласа сигналу на вході при нульових початкових умовах:

(1)

Оскільки відомі ПФ всіх елементів, що входять в структурну схему (рис.1), то застосовуючи апарат структурних перетворень, що дозволяє знаходити ПФ замкнутих систем, заданих структурними схемами, отримаємо ПФ розімкнутої та замкнутої САУ, зображеної на рис.1:

(1)

У формулу (1) підставлені чисельні значення, взяті з таблиці 1.

(2)

Отримано ПФ замкнутої системи (2).

2.2 Знаходження коренів характеристичного рівняння, використовуючи МПІ

Для того, щоб лінійна стаціонарна система була стійкою, всі корені її характеристичного рівняння (полюси ПФ) повинні розташовуватися в лівій половині s-площині.

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

2.2.1 Вид характеристичного рівняння

Запишемо характеристичне рівняння знайденої ПФ (формула 2):

2.2.2 Метод січних.

Проведемо локалізацію коренів:

Побудуємо графік функції на інтервалі :

Рис.2. Графік характеристичного полінома (3) на інтервалі

Рівняння має 1 дійсний корінь і 2 уявних.

Рівняння вирішується методом січних (4):

(4)

Візьмемо початкове наближення і для знаходження дійсного кореня.

S =- 8.210097

Далі отримаємо значення комплексних коренів:

Підставимо в (5)

Отримуємо коріння характеристичного рівняння:

Висновок: 2 полюси передавальної функції знаходяться в правій півплощині. Система нестійка.

2.2.3 Рух дійсного кореня полінома в s-площині

Побудуємо графік руху кореня в залежності від номера ітерації:

Рис.3. Графік руху кореня в залежності від номера ітерації

2.3 Аналітичні вирази для АЧХ, ФЧХ, АФЧХ

Графік АЧХ:

Функції, зумовлені залежностями (6) і (7), називаються відповідно амплітудно-частотної (АЧХ) і фазочастотного (ФЧХ) характеристиками.

Частотні характеристики визначаються наступними показниками:

показник коливальності - Характеризує схильність системи до коливань: чим вище , Тим менш якісна система (як правило в реальних системах );

резонансна частота - Частота, при якій АЧХ має максимум (на цій частоті гармонійні коливання мають найбільше посилення);

смуга пропускання системи - інтервал від до , При якому виконується умова ;

частота зрізу - Частота, при якій АЧХ системи приймає значення, рівне , Тобто ;

Частота зрізу побічно характеризує тривалість перехідного процесу; справедливе співвідношення . Таким чином можна зробити висновок: чим ширше смуга пропускання, тим система є більш швидкодіючою.

2.4. Годограф АФЧХ і графіки АЧХ і ФЧХ із зазначенням частот

Рис.4 Графік АЧХ заданої САУ

Рис.5 Графік ФЧХ заданої САУ

Рис.6 Графік АФЧХ заданої САУ

2.5 Диференціальне рівняння заданої САУ

Отримаємо ДУ заданої САУ:

2.6 Нормальна форма Коші, отриманого ДУ 3-го порядку

Так як ДУ заданої САУ має вищий порядок, то його необхідно звести до системи рівнянь, кожне з яких повинен мати перший порядок, тобто має місце нормальна форма Коші:

. (9)

Так як ДУ заданої САУ має укорочену праву частину, то запишемо нормальну форму Коші в наступному вигляді:

. (10)

Наведемо рівняння (12) до нормальної формі Коші:

(11)

або

де

2.7 Аналітичне рішення ДУ

Нехай задано зображення виходу

або .

Тоді використовуючи другу теорему розкладання Лапласа отримаємо наступне аналітичний вираз для вихідного сигналу:

реакція системи на одиничне поетапне вплив ( ) (12):

(12)

2.8 Рішення ДУ чисельним методом (метод Рунге-Кутта 5-го порядку і метод Адамса неявний 4-го порядку)

У неявних методах використовується інформація про можливе майбутньому значенні рішення в точці п +1. Це дещо підвищує точність одержуваних результатів у порівнянні з явними методами.

Для організації обчислювального процесу по інтерполяційної формулою Адамса, що має точність рішення (13):

необхідно заготовити початкові значення , Використовуючи метод Рунге-Кутта 5-її порядку.

Наведені коефіцієнти:

Проведемо дослідження рішення ДУ в залежності від кроку:

Графіки вихідного сигналу, отриманого в аналітичному вигляді, вихідного сигналу, отриманого рішенням ДК та похибки рішення при кроці h = 0.1 і h = 0.01, h = 0.001.

Рис.7. Графіки вихідного сигналу , Отриманого в аналітичному вигляді, вихідного сигналу , Отриманого чисельним рішенням ДК та похибки рішення при кроці

Рис.8. Графіки вихідного сигналу , Отриманого в аналітичному вигляді, вихідного сигналу , Отриманого чисельним рішенням ДК та похибки рішення при кроці

Рис.9. Графіки вихідного сигналу , Отриманого в аналітичному вигляді, вихідного сигналу , Отриманого чисельним рішенням ДК та похибки рішення при кроці

2.9 Аналізу заданої системи з використанням спектрального методу (базис: Чебишева 2 роди)

Спектральна форма представлення сигналів і тимчасових динамічних характеристик систем та об'єктів заснована на їх розкладанні в заданій системі ортогональних функцій

Якщо певний сигнал належить простору , Тобто для нього справедливо положення

То він може бути представлений у вигляді ряду Фур'є:

(14)

Якщо ввести вектори

то ряд (14) можна представити таким чином

(15)

Сукупність коефіцієнтів Фур'є розкладання сигналу в ряд (14) називається спектральною характеристикою цього сигналу.

Коефіцієнти Фур'є визначаються за формулою

(16)

Суттєвим і визначальним відзнакою спектрального опису дискретних сигналів від спектрального опису безперервних сигналів на кінцевих інтервалах є можливість їх точного уявлення у вигляді рядів Фур'є з кінцевим числом членів. Значить, якщо дискретний сигнал, а даний сигнал має місце на вході ЕОМ після його аналого-цифрового перетворення (АЦП), заданий на кінцевій множині точок, наприклад , У вигляді деякої числової послідовності , То його розкладання по заданій системі ортогональних функцій

визначається співвідношенням

(17)

Система - Це система ортогональних, нормованих функцій, що задовольняють умові

Коефіцієнти Фур'є визначаються за формулою

(18)

Далі вводимо поліноми Чебишева 2-го роду (19):

(19)

2.9.1 Алгоритм побудови спектральної характеристики (СХ)

1. Вихідні рівняння (20):

(20)

Обчислимо ядра і (21):

(21)

3. Розкладемо в ряди Фур'є по заданому базису (22):

(22)

4. Отримаємо значення Сх з наведених нижче перетворень (23):

(23)

5. Знайдемо матрицю А:

6. Отримані значення ядер:

7. Дія:

8. Значення вектора Cх:

9. Матриця А:

А =

Рис.10 Перехідна функція, побудована спектральним методом

Рис.11 Графік вихідного сигналу, отриманого аналітично, сигналу, отриманого спектральним методом і помилки.

3. СИНТЕЗ

Вихідні дані: структурна схема заданої системи зображена на рис. 12.

Введемо в систему послідовне коригуючий пристрій. Як регулятор виберемо ПІД-регулятор.

Його передатна функція має вигляд:

(24)

Рис.12: Структурна схема заданої САУ з коригувальним пристроєм у прямій ланцюга.

3.1 Передавальна функція замкнутої ланцюга скоригованого САУ

Знайдемо передавальну функцію розімкнутого ланцюга, якщо відома передатна функція об'єкта (25):

(25)

Передавальна функція замкнутої системи буде мати вигляд (26):

(26)

Для вирішення задачі синтезу необхідно знайти параметра регулятора , Структура якого заданна (формула 31), при яких реальний вихідний сигнал , Що є реакцією на одиничне поетапне вплив, буде близький до заданого еталонному сигналу .

В якості еталонного вихідного сигналу виберемо наступний сигнал:

, (27)

де параметр знаходиться за такою залежністю:

. (28)

3.2 Функціонал якості, який підлягає подальшій мінімізації

Критерієм близькості виберемо метрику простору .

Тоді цільова функція, що підлягає мінімізації за параметрами регулятора буде мати наступний вигляд:

(29)

3.2.1 Пошук мінімуму функції методом Фібоначчі

Якщо початковий інтервал має довжину , То провівши обчислень функції, можна зменшити початковий інтервал невизначеності в раз за такою формулою:

(30)

в порівнянні з його початковою довжиною (нехтуючи ).

Якщо визначити послідовність чисел Фібоначчі наступним чином: для то можна знайти положення першої точки, яка знаходиться на відстані від одного з кінців початкового інтервалу, причому не важливо, від якого кінця, оскільки друга точка поміщається згідно з правилом симетрії на відстані від кінця інтервалу:

. (31)

Після того, як знайдено положення першої точки, числа Фібоначчі більше не потрібні. Використовуване значення може визначаться з практичних міркувань. Воно повинно бути менше , Інакше будуть мати місце зайві обчислення значень функції .

Таким чином, пошук методом Фібоначчі є ітераційної процедури.

У процесі пошуку інтервалу з точкою , Вже лежить в цьому інтервалі, наступна точка завжди вибирається такий, що .

Позначимо і , Тоді можна розглянути чотири випадки організації обчислювального процесу:

1. : Новий інтервал .

2. : Новий інтервал .

3. : Новий інтервал .

4. : Новий інтервал .

Закінчувати обчислювальний процес можна двома способами. Або виконати намічені раніше обчислень, або, якщо в процесі обчислень інтервал невизначеності стане менше заданої величини.

3.2.2 Метод Хука-Дживса

У даному методі пошук складається з послідовності кроків досліджує пошуку навколо базисної точки, за яким, у разі успіху, проводиться пошук за зразком.

Процедура пошуку наступна.

Вибрати початкову базисну точку і крок довжиною для кожної з змінних , .

Обчислити у базисній точці з метою отримання відомостей про локальному поведінці функції . Ці відомості будуть використовуватися для знаходження відповідного напрямку пошуку за зразком, за допомогою якого можна сподіватися досягти більшого зменшення значення функції.

При пошуку за зразком використовується інформація, отримана в процесі дослідження, і мінімізація функції завершується пошуком в напрямку, заданому зразком.

Завершити цей процес, коли довжина кроку (довжини кроків) буде зменшуватися до заданого малого значення.

3.3 Диференціальне рівняння скоригованої системи

Для мінімізації цільової функції (37) необхідно реалізувати обчислення реального вихідного сигналу в кожен окремий момент часу. Крім цього, необхідно реалізувати ітераційний процес і реалізувати алгоритм обчислення параметра

Для обчислення перейдемо від передавальної функції замкненого кола до диференціальних рівнянь, використовуючи властивості перетворень Лапласа. Воно буде мати наступний вигляд:

(32)

Запишемо ДУ (32) в іншому вигляді:

(33)

3.4 Нормальна форма Коші, отриманого ДУ скоригованої системи

Для вирішення ДУ (33) за допомогою чисельного методу розв'язання диференціальних рівнянь, необхідно знизити його порядок, шляхом переходу від даного ДУ до нормальної форми Коші

Нормальна форма Коші для ДУ (33) буде мати наступний вигляд:

де коефіцієнти розраховуються за такими формулами:

Тоді ДУ (33) можна записати в наступному вигляді

, (34)

де

Рис. 13. Графіки вихідних сигналів скоригованого (зелена лінія) і не скоригованого (синя лінія) САУ.

Отримані параметри регулятора:

Кп = 1.0547895

Кд = 0.0550905

Ки = 0.9452075

5. Висновки

Чисельні методи рішення диференціальних рівнянь використовуються в тих випадках, коли не вдається знайти їх рішення в аналітичному вигляді. Перш за все, це відноситься до лінійних диференціальних рівнянь зі змінними коефіцієнтами та нелінійним диференціальним рівнянням, відповідно описує динаміку лінійних нестаціонарних і нелінійних систем управління.

Сутність чисельних методів полягає в тому, що рішення ДУ будується тільки для дискретних значень аргументу.

Всі чисельні рішення ДУ діляться на дві групи: однокрокові і багатокрокові. У однокрокових методах використовується інформація про поведінку рішення в попередній точці. У багатокрокових про поведінку рішення в декількох попередніх точках.

Чисельні рішення ДУ можна розділити на дві групи: явні і неявні. В явних методах, на відміну від неявних, використовується явна залежність значення функції в поточній крапці від значень функції в попередніх точках. Перевагою таких методів є відносна простота обчислення значення функції на кожному кроці, проте, збіжність даних методів визначається кроком інтегрування .

У відношенні чисельних методів оптимізації слід зазначити наступне. Всі чисельні методи мінімізації діляться на прямі і градієнтні методи. У прямих методах використовується тільки значення функції в конкретних точках, а в градієнтних - інформація про перших і других похідних функції. Також методи мінімізації можна розділити на методи мінімізації функції однієї змінної і методи, що дозволяють мінімізувати функції багатьох змінних. При мінімізації необхідно враховувати наявність обмежень на параметри вихідної функції.

6. Література

Вержбицький. Чисельні методи.

Конспект лекцій з курсу "Системи аналітичних обчислень" за I і II семестр.

7. Додаток 1 (Лістинг скриптів для знаходження коренів полінома)

function secush

clc

e = 10. ^ -5;

x =- 8.1;

xm1 =- 8

Asm1 = 8.6159999

i = 0;

As = 0.252 * (x. ^ 3) +1.41 * (x. ^ 2) +14.2 * x +161;

x1 = x-(As. * (xm1-x ))./( Asm1-As);

Asm1 = As;

As = 0.252 * (x1. ^ 3) +1.41 * (x1. ^ 2) +14.2 * x1 +161;

i = i +1;

while abs (x1-x)> e

xm1 = x;

x = x1;

x1 = x-(As. * (xm1-x ))./( Asm1-As);

Asm1 = As;

As = 0.252 * (x1. ^ 3) +1.41 * (x1. ^ 2) +14.2 * x1 +161;

i = i +1;

A (i) = x;

end

hold on

for n = 1: i

plot (n, A (n), 'b-o')

end

grid on

xlabel ('iteraciya')

ylabel ('roots')

disp ('відповідь');

disp (x);

8. Додаток 2 (Лістинг скриптів для вирішення ДУ)

function Difer

clc

T = 4;

a0 = 638.89;

a1 = 56.35;

a2 = 5.60;

b0 = 595.24;

h = 0.0005;

A_X (1,1:3) = [0 0 0];

A = [0 1 0;

0 0 1;

a0 a1 a2];

B = [0 0 b0] ';

k = 0;

t = 0;

while (t <(Th))

if (t <= 3 * h)

K1 = A * (A_X (k +1 ,:))';

K2 = A * (A_X (k +1 ,:))'+ 1 / 3 * K1;

K3 = A * (A_X (k +1 ,:))'+ 1 / 6 * K1 +1 / 6 * K2;

K4 = A * (A_X (k +1 ,:))'+ 1 / 8 * K1 +3 / 8 * K2;

K5 = A * (A_X (k +1 ,:))'+ 1 / 2 * K1-3 / 2 * K3 +2 * K4;

A_X (k +2 ,:)=( A_X (k +1 ,:))+ h / 6 * (K1 '+4 * K4' + K5 ');

else

h1 = h;

t = t + h1;

H = (eye (length (A_X (1 ,:)))-( 9 * h1/24) * A);

G = (eye (length (A_X (1 ,:)))+ 19 * h1/24 * A) * (A_X (k +1 ,:))'+ h1/24 * A * (-5 * (A_X ( k ,:))'+( A_X (k-1 ,:))')

+ H1/24 * B * (9 * 1 +19 * 1-5 * 1);

A_X (k +2 ,:)=( inv (H) * G) ';

end

Otr (k +1) = t;

k = k +1;

h =- 0.43496

end

plot (Otr, A_X (1: k, 1), 'b-');

grid on

9. Додаток 4 (Лістинг скриптів для спектрального аналізу)

spectr.m

syms t T;

Kx = (638.89 / 2) * (tT). # 2-56.35 * (1. / 2) * (-2 * (tT)) +5.6;

Ky = (595.24 / 2) * (tT). # 2;

F2 = 2 * t;

L (2) = F2;

F3 = 4 * t. ^ 2-1;

L (3) = F3;

F4 = 8 * t. ^ 3-4 * t;

L (4) = F4;

F5 = 16 * t. ^ 4-12 * t. ^ 2 +1;

L (5) = F5;

F6 = 32 * t. ^ 5-32 * t. ^ 3 +6 * t;

L (6) = F6;

F7 = 64 * t. ^ 6-80 * t. ^ 4 +24 * t. ^ 2-1;

L (7) = F7;

F8 = 128 * t. ^ 7-192 * t. ^ 5 +80 * t. ^ 3-8 * t;

L (8) = F8;

F9 = 256 * t. ^ 8-448 * t. ^ 6 +240 * t. ^ 4-40 * t. ^ 2 +1;

L (9) = F9;

F10 = 512 * t. ^ 9-1024 * t. ^ 7 +672 * t. ^ 5-160 * t. ^ 3 +10 * t;

L (10) = F10;

F1 = 1;

L (1) = F1;

F2 = 2 * T;

L1 (2) = F2;

F3 = 4 * T. ^ 2-1;

L1 (3) = F3;

F4 = 8 * T. ^ 3-4 * T;

L1 (4) = F4;

F5 = 16 * T. ^ 4-12 * T. ^ 2 +1;

L1 (5) = F5;

F6 = 32 * T. ^ 5-32 * T. ^ 3 +6 * T;

L1 (6) = F6;

F7 = 64 * T. ^ 6-80 * T. ^ 4 +24 * T. ^ 2-1;

L1 (7) = F7;

F8 = 128 * T. ^ 7-192 * T. ^ 5 +80 * T. ^ 3-8 * T;

F9 = 256 * T. ^ 8-448 * T. ^ 6 +240 * T. ^ 4-40 * T. ^ 2 +1;

L1 (9) = F9;

F10 = 512 * T. ^ 9-1024 * T. ^ 7 +672 * T. ^ 5-160 * T. ^ 3 +10 * T;

L1 (10) = F10;

F1 = 1;

L1 (1) = F1;

G = L '* L1;

In = Kx * G;

r = int (In, T, 0, t);

Cx = int (r, t, 0,1.5);

In = Ky .* G;

r = int (In, T, 0, t);

Cy = int (r, t, 0,1.5);

A = ((Cx + eye (10 )).^- 1) * Cy;

Cy = int (L, t, 0,1.5);

Cx = A * Cy '

Postr.m

function H = fun (t)

Cx = [3.7672; 1.3134; 0.5181; 0.2065; 0.0819; 0.0323; 0.0127; 0.0491; 0.0189; 0.0723];

F2 = 2 * t;

L (2) = F2;

F3 = 4 * t. ^ 2-1;

L (3) = F3;

F4 = 8 * t. ^ 3-4 * t;

L (4) = F4;

F5 = 16 * t. ^ 4-12 * t. ^ 2 +1;

L (5) = F5;

F6 = 32 * t. ^ 5-32 * t. ^ 3 +6 * t;

L (6) = F6;

F7 = 64 * t. ^ 6-80 * t. ^ 4 +24 * t. ^ 2-1;

L (7) = F7;

F8 = 128 * t. ^ 7-192 * t. ^ 5 +80 * t. ^ 3-8 * t;

L (8) = F8;

F9 = 256 * t. ^ 8-448 * t. ^ 6 +240 * t. ^ 4-40 * t. ^ 2 +1;

L (9) = F9;

F10 = 512 * t. ^ 9-1024 * t. ^ 7 +672 * t. ^ 5-160 * t. ^ 3 +10 * t;

L (10) = F10;

F1 = 1;

L (1) = F1;

H = (Cx '* L');

t = [0:0.01:5];

plot (t, H)

10. Додаток 5 (Лістинг скриптів для оптимізації)

jivs.m

clear

clc

a = 0;

b = 5;

h = 0.1;

Kp1 = 2; Kd1 = 1; Ki1 = 0;

J1 = int2 (Kp1, Kd1, Ki1);

while h> 0.0000001

Kp2 = Kp1 + h;

J2 = int2 (Kp2, Kd1, Ki1);

if J2> J1

Kp2 = Kp1-h;

J2 = int2 (Kp2, Kd1, Ki1);

if J2> J1

Kp2 = Kp1;

end

Kd2 = Kd1 + h;

J2 = int2 (Kp2, Kd2, Ki1);

if J2> J1

Kd2 = Kd1-h;

J2 = int2 (Kp2, Kd2, Ki1);

if J2> J1

Kd2 = Kd1;

end

Ki2 = Ki1 + h;

J2 = int2 (Kp2, Kd2, Ki2, h);

if J2> J1

Ki2 = Ki1-h;

J2 = int2 (Kp2, Kd2, Ki2, h);

if J2> J1

Ki2 = Ki1;

end

h = fibon (a, b, h);

while J2 <J1

Kp = Kp1 +2 * (Kp2-Kp1); Kd = Kd1 +2 * (Kd2-Kd1); Ki = Ki1 +2 * (Ki2-Ki1);

J1 = J2;

J2 = int2 (Kp, Kd, Ki, h);

Kp1 = Kp2; Kp2 = Kp; Kd1 = Kd2; Kd2 = Kd; Ki1 = Ki2; Ki2 = Ki;

end

disp (Kp)

disp (Kd)

disp (Ki)

int2.m

function J = int2 (Kp, Kd, Ki, h)

clc

T = 4;

A = [0 1 0 0; 0 0 1 0; 0 0 0 1; -595.23809523 * Ki -43.6507936507-595.23809523 * Kp -56.34920634920635-595.23809523 * Kd -5.59523809523809];

B = [0; 595.23809 * Kd; 595.23809 * Kp-3330.498866 * Kd; 595.23809 * Ki-33540.615-354308.277 * (Kd) ^ 2-3330.498 * Kp-18634.934 * Kd];

k = 0;

t = 0;

while (t <(Th))

if (t <= 3 * h)

K1 = A * (A_X (k +1 ,:))';

K2 = A * (A_X (k +1 ,:))'+ 1 / 3 * K1;

K3 = A * (A_X (k +1 ,:))'+ 1 / 6 * K1 +1 / 6 * K2;

K4 = A * (A_X (k +1 ,:))'+ 1 / 8 * K1 +3 / 8 * K2;

K5 = A * (A_X (k +1 ,:))'+ 1 / 2 * K1-3 / 2 * K3 +2 * K4;

A_X (k +2 ,:)=( A_X (k +1 ,:))+ h / 6 * (K1 '+4 * K4' + K5 ');

else

h1 = h;

t = t + h1;

H = (eye (length (A_X (1 ,:)))-( 9 * h1/24) * A);

G = (eye (length (A_X (1 ,:)))+ 19 * h1/24 * A) * (A_X (k +1 ,:))'+ h1/24 * A * (-5 * (A_X ( k ,:))'+( A_X (k-1 ,:))')

+ H1/24 * B * (9 * 1 +19 * 1-5 * 1);

A_X (k +2 ,:)=( inv (H) * G) ';

end

Otr (k +1) = t;

k = k +1;

end

grid on

fibon.m

function h = fibon (a, b, h)

F (1) = 1; F (2) = 1; n = 100;

for i = [1:0.1: n-2]

F (i +2) = F (i +1) + F (i);

end

j = 0;

x1 = a; x3 = b;

L1 = x3-x1;

L2 = (F (n-1) / F (n)) * L1 + ((-1) ^ n) / F (n) * eps;

x2 = x3-L2;

x4 = x1 + x3-x2;

while (abs (x3-x1)> eps)

F2 = x2;

F4 = x4;

if ((x2 <x4) & & (norm (F2) <norm (F4)))

x1 = x1; x3 = x4;

x4 = x1 + x3-x2;

elseif ((x2> x4) & & (norm (F2) <norm (F4)))

x1 = x4; x3 = x3;

x4 = x1 + x3-x2;

elseif ((x2 <x4) & & (norm (F2)> norm (F4)))

x1 = x2; x3 = x3;

x2 = x1 + x3-x4;

elseif ((x2> x4) & & (norm (F2)> norm (F4)))

x1 = x1; x3 = x2;

x2 = x1 + x3-x4;

end

j = j +1;

la = x1 + (x3-x1) / 2;

end

l = la;