При організації чисельних експериментів у SIMULINK може виникнути необхідність багаторазово запускати ту ж саму модель з різними параметрами і вхідними сигналами. Це можна забезпечити, керуючи процесом моделювання з командного вікна. Насамперед, необхідно відкрити Simulink-модель вручну чи за допомогою команди open. Так, команда open('vdp'). відкриває модель vdp.mdl (схема моделювання рівняння Вандерполя). Щоб запустити процес моделювання з командного рядка використовується команда sim, наприклад sim('vdp'). Як параметри цій команді можна передати час моделювання, вхідний сигнал (при наявності вхідних портів у моделі), початкові умови й опції решателя. У цьому випадку синтаксис її виклику буде мати вид
y=sim('model',timespan,options,ut),
де timespan – інтервал моделювання, ut – масив, що складається з отсчетов часу і значень вхідного сигналу.
Приклад 3.Створимо модель, показану на рис., і збережемо її під ім'ям simple_mo-del. Знайдемо її реакцію на синусоїду на інтервалі 30 с. Опцій моделювання використовувати не будемо – як аргумент візьмемо порожній масив.
Після виконання цих команд можна відкрити блок Scopeу схемі і побачити результати моделювання. Помістити результати моделювання в робочий простір і потім подивитися графік за допомогою команди plotможна, використовуючи виклик simз вихідними параметрами.При виклику команди simз опціями решателя структура опцій задається командою simset. Для приклада змінимо тип решателя на ode23:
З повним списком опцій моделювання можна ознайомитися, використовуючи команду simget, наприклад o=simget ('simple_model').Крім завдання опцій моделювання, існує можливість програмно змінювати параметри блоків.
Для їхнього читання і запису використовують команди get_paramі set_param. Оскільки назви параметрів блоків не завжди збігаються з текстами пояснень у вікні параметрів, буває корисна конструкція get_param(gcb, 'objectparameters'). Мишею у вікні моделі вибирається потрібний блок. Після цього команда gcbповертає ім'я поточного блоку.
Так, якщо буде виділений блок Gain, команда gcbповерне рядок 'simple_model/ Ga-in'. За допомогою тієї ж функції get_paramможна одержати будь-як конкретний параметр, наприклад:
Зверніть увагу на те, що тип параметра – рядок. Це значить,що якщо, наприклад, необхідно змінити коефіцієнт підсилення з 3 на 4 за допомогою функції set_param, то прийдеться набрати не set_param('simple_model/Gain', 'Gain',4), а set_param('simple_model/Gain', 'Gain', num2str(4)).Використання цих команд дозволяє, зокрема, організувати циклічну зміну параметрів схеми моделювання (коефіцієнтів підсилення, початкових умов) за допомогою оператора MATLAB for.
Приклад 4(задача про оптимальні початкові умови). Рух маятника описується нелінійним диференціальним рівнянням . Його початкові умови задовольняють умові (лежать на колі радіуса R у фазовій площині). Потрібно знайти значення , при яких енергія вихідного сигналу буде максимальної або мінімальною: .
Рішення. Складемо схему моделювання в SIMULINK і назвемо її sys.mdl (рис. ). Розглянемо два способи пошуку оптимальних початкових умов.
Спосіб 1 (одновимірний лінійний пошук).Будемо запускати схему з різних початкових умов, обчислюючи щораз енергію Е и запам'ятовуючи результат, а потім виберемо екстремальні значення.
Рис. Sys.mdl
Початкові умови будемо формувати в MATLAB за формулами Програму для керування процесом моделювання оформимо у вигляді окремого m_файлу optic.m.
Program optic %Optimal Initial Condition for sys.mdl
open('sys'); % виклик моделі sys.mdl
E=zeros(1,200); R=1; % масив для запису енергій
for і=1:200
x0=R*cosd(і); y0=R*sind(і); % чергові початкові умови
set_param('sys/х','initial',num2str(x0)); % установка початкових умов
set_param('sys/у','initial',num2str(y0));
sim('sys',[0,40],[],[]); % запуск моделювання
E(i)=trapz(tout,yout.^2); % обчислення енергії
disp(i); % висновок поточного i на екран
end
plot(E,'LineWidth',2);grid % графік залежності енергії
title('ENERGY vs. IC angle'), % від кута на одиничній окружності
Результат моделювання для R = 1 показаний на рис.
З графіка видно, що максимальна енергія вихідного сигналу приблизно дорівнює 2,4.Більш точні дані одержуємо, додаючи в програму наступні команди:
[EM,iM=max(E);
x0M=cosd(iM);
y0M=sind(iM);
disp('im,x0m y0m Em='),[im,x0m y0m Em],
[Em,im]=min(E);
x0m=cosd(im);
y0m=sind(im);
disp('iM,x0M y0M EM='),[iM,x0M y0M EM],
Результат їхнього виконання:
iM,x0M y0M EM = 47 0.6820 0.7314 2.4118
im,x0m y0m Em = 136 _0.7193 0.6947 1.2637
Таким чином, максимальна енергія вихідного сигналу, рівна Е = 2,412, досягається при j= 47про, тобто при початкових умовах (маятник відхиляємо вправо і штовхаємо від положення рівноваги). Мінімальна енергія вихідного сигналу,рівна Е = 1,264, досягається при j = 136про, тобто при початкових умовах маятник відхиляємо вліво і штовхаємо до положення рівноваги). Відповідні графіки приведені на рис. Б (вгоры).
Змінимо радіус R окружності початкових умов, узявши його близьким до критичного значення, при якому кут у відхилення маятника наближається до p. У цьому випадку маятник починає «зависати» у верхнім положенні, що приводить до різкого збільшення енергії вихідного сигналу ( ). На рис. а (внизу) показані графіки енергії Е в залежності від кута _ і максимального і мінімального вихідного сигналу в у залежності від часу (рис., б) для R = 2,29, там же приведені чисельні значення екстремальних параметрів.
Спосіб 2 (пошук екстремума за допомогою функції fmincon- рішення задач на умовний і безумовний екстремум. У даному випадку треба знайти екстремуми функції Е (енергії вихідного сигналу), що залежить від двох аргументів (початкових умов маятника) при наявносты обмеження (сума квадратів початкових умов фіксована):
Для рішення таких задач призначена команда fmincon.Для її застосування потрібно створити дві допоміжні функції, що містять опис критерію, що мінімізується, й обмежень. Оформимо ці функції у вигляді m-файлів з іменами
energyі nlcon:
Запускаємо процес оптимізації, набираючи в робочому вікні команду fmincon:
По закінченні ітераційного процесу одержуємо повідомлення про його успішне завершення і значення вихідних параметрів (вектор початкових умов і мінімальну енергію для R = 1). Щоб вирішити задачу про максимальну енергію, потрібно змінити знак перед командою trapz в останньому рядку функції energyі повторити обчислення. У результаті одержуємо
Знайдені значення близькі до отриманих першим способом, але більш точні. Виконуючи аналогічні обчислення для R = 2.29, також одержуємо уточнені значення оптимальних початкових умов:
X = 1.7201 2.2900 E =-51.8823 X = 1.2100 2.2900 E = 6.5955
1.5118 - _1.9442
Звернемо увагу на два нюанси в оформленні функції energy.Якщо в другому рядку забрати знак коментарю (%), то виклик моделі sys.mdl і перемальовування схеми моделювання буде відбуватися на кожнім кроці ітераційного процесу, що сильно збільшить час
Рі-шення. Простіше відкрити цю модель вручну до початку рішення, тоді команда open('sys') не буде потрібно. Другий нюанс стосується синтаксису команди num2str(x0,15). Її другий аргумент задає кількість розрядів числа x0, преутворених у рядок. Якщо його опустити,то за замовчуванням воно може виявитися рівним 4, що цілком спотворить результати оптимізації.
З подібними тонкостями нерідко приходиться зустрічатися при моделюванні в MATLAB і інших пакетах. Вони вимагають від користувача гарного розуміння сутності розв'язуваної задачі і застосовуваного математичного апарата з одного боку, і знання мови програмування і розуміння процесу комп'ютерної реалізації використовуваного алгоритму – з іншого
ІІ
РОБОТА З РОЗКИДАНИМИ ДАНИМИ
У MATLAB входять функції для роботи з розкиданими даними, тобто наприклад, у двовимірному випадку це дані zi, задані в крапках (xi, yi)i=1,2,...,k. Ці функції дозволяють вирішувати наступні задачі.
· Інтерполювати розкидані дані різними способами. Для даних, визначених у крапках на площині, можливо кусочно-постійне наближення по найближчому сусіді, кусочно-лінійне наближення і більш гладке наближення з неперервними похідними (функція griddata). Дані, визначені в крапках тривимірного і n-мірного простору можуть бути наближені по найближчому сусіді і за допомогою кусочно-лінійних функцій (функції griddata3 і griddatan). Значення интерполянта повертається на ортогональній сітці, що накриває крапки, у яких задані розкидані дані.
· Будувати опуклу оболонку множини крапок на площині, у тривимірному і n-мірному просторі (функції convhull, convhulln).
· Реалізація опуклої оболонки множини крапок і її візуалізація
· Опукла оболонка conv(P) множини крапок, що належать n-мірному простору Rn, визначається як найменша опукла множина, що містить усі ці крапки. Опукла множина, породжена крапками (pi)i=1,2,...,k, визначається як множина крапок:
· p = a1p1 + a2p2 + ... + akpk
· де всі ai >/= 0 і a1 + a2 + ... + ak = 1 (тут крапки ототожнюються з векторами, що виходять з початку координат).
· У двовимірному випадку є проста інтуїтивно зрозуміла аналогія. Нехай на площині лежать крапки (убиті цвяхи), ми натягаємо навколо них мотузку.
·
· Те, що вийде, і буде опуклою оболонкою (на правому малюнку її границя зображена червоним кольором). Усі крапки з множини знаходяться або усередині опуклої оболонки, або на її границі.
· Опукла оболонка є перетинанням кінцевого числа півпросторів (півпростір - частина простору, що лежить по одну сторону від гіперплощини). У двовимірному випадку гіперплощини це лінії, а в тривимірному випадку - площини.
·
Для побудови опуклої оболонки в MATLAB служить функція convhull. Її вхідними аргументами є вектори з координатами абсцис і ординат крапок, а вихідними аргументами є номери крапок, що утворять вершини опуклої оболонки. У наступному прикладі генерується набір з 100 крапок з координатами, що змінюються випадковим образом від 0 до 1, потім за допомогою функції convhull будується опукла оболонка. Номера крапок, що є вершинами опуклої оболонки, заносяться в масив kvert. Потім вершини опуклої оболонки виділяються круглими маркерами і малюється сама опукла оболонка
Якщо функцію convhull викликати з другим додатковим вихідним аргументом, то в нього запишеться площа опуклої оболонки
>> [kvert, s] = convhull(x, y);
>> s
s = 0.9164
Для побудови опуклої оболонки крапок із тривимірного простору і просторів більш високих вимірностей використовується функція convhulln. Координати крапок задаються в рядках двовимірного масиву, що вказується як вхідний аргумент функції convhulln. Вона повертає двовимірний масив, кожен рядок якого містить номера крапок з вихідної множини, що утворять грань опуклої оболонки (у двовимірному випадку грань - відрізок, що з'єднує сусідні вершини опуклої оболонки).
Як приклад приведемо побудову опуклої оболонки крапоок, що випадковим образом розподілені в кубі [0,1]x[0,1]x[0,1]. Координати крапок генеруються за допомогою функції rand і записуються у вектори-стовпці. З цих стовпців створюється матриця, що вказується як вхідний аргумент функції convhulln. Функція convhulln повертає матрицю з трьох стовпців, кожний з який містить номера крапок вихідної множини, що утворять деяку грань опуклої оболонки. Для візуалізації опуклої оболонки застосовується функція trimesh, що будує полігональний об'єкт (patch) і повертає покажчик на нього. Вид цього полігонального об'єкта можна змінювати за своїм розсудом, включаючи прозорість і колір граней, маркери у вершинах і т.д. Коли опукла оболонка знайдена, та уся вихідна множина крапок розбивається на дві підмножини - крапки, що є вершинами опуклої оболонки, і внутрішні крапки.
У прикладі, що нижче приводиться, поруч із крапками цих типів розміщаються підписи різного кольору з номерами крапок. Для цього використовується функція text. При цьому виникає питання, як визначити номера крапок, що потрапили усередину опуклої оболонки вихідної множини, і тих, котрі належать вершинам опуклої оболонки. Як уже говорилося, номера вершин записані у вихідному аргументі функції convhulln (двовимірному масиві з трьома стовпцями). Однак, вони зустрічаються в різних рядках цього масиву і не по одному разі. Тому для виділення номерів крапок, що потрапили у вершини опуклої оболонки, можна робити так: послідовно об'єднати вміст стовпців матриці з номерами вершин, скориставшись функцією union. Оскільки дана функція робить об'єднання в змісті об'єднання множин, то повтори елементів виключаються. Аналогічно, для перебування номерів крапок, що потрапили усередину опуклої оболонки, можна застосувати функцію setdiff, що знаходить різницю множин. У нашому випадку з множини номерів усіх крапок варто відняти номера крапок, що утворять вершини опуклої оболонки.
k=20; задаємо двадцять крапок
rand('state',0); генератору випадкових чисел задаємо початковий стан
% генеруємо координати крапок множини і записуємо їх у вектори-стовпці x, y, і z
x=rand(k,1);
y=rand(k,1);
z=rand(k,1);
% утворимо двовимірний масив з координатами крапок множини
X=[x y z];
% знаходимо номера крапок, що утворять трикутної грані опуклої оболонки множини
kvert=convhulln(X)
% визуализируем опуклу оболонку за допомогою функції trisurf
figure('Color','w')
h=trisurf(kvert,X(:,1),X(:,2),X(:,3))
hold on
% установлюємо прозорість граней, червоні ребра
% і круглі червоні маркери у вершинах опуклої оболонки
set(h,'FaceAlpha',0.1,'FaceColor','g',...
'EdgeColor','r','Marker','о','MarkerSize',10)
% знаходимо номера крапок, що є вершинами опуклої оболонки, і записуємо їх у kv
У результаті виконання цих команд на осі графічного вікна виводиться опукла оболонка, її пронумеровані вершини і пронумеровані крапки, що потрапили усередину опуклої оболонки, відповідно до їхніх номерів у вихідній множини:
Замість функції trisurf можна використовувати функцію trimesh, що на відміну від trimesh не заливає грані кольорами.
При виклику функції convhulln(X) із другим додатковим вихідним аргументом, у ньому повертається обсяг опуклої оболонки
>> [kvert,V]=convhulln(X);
>> V
V =
0.2582
У MATLAB для побудови опуклої оболонки використовується функція qconvex пакета Qhull (http://www.qhull.org).Ці функція підтримує ряд опцій для завдання деяких параметрів алгоритму побудови опуклої оболонки, який можна задавати й у функціях convhull, convhulln.
Алгоритми побудови опуклої оболонки, у тому числі й алгоритм Quickhull, описані, наприклад, на http://algolist.manual.ru/maths/geom/convhull/. Там же приведено кілька посилань на статті з описами алгоритмів побудови опуклої оболонки множини крапок.
Для візуалізації опуклої оболонки крапок тривимірного простору не обов'язково використовувати функцію trisurf. Замість її можна самому сконструювати полігональний об'єкт із заданими способами заливання кольором ребер і граней і їхньою прозорістю.
Одним зі способів завдання полігонального об'єкта є вказівка його вершин і матриці, що містить номери вершин кожної грані, тобто саме тієї матриці, що повертається функцією convhulln. При цьому, зокрема, можна вказати колір кожної вершини і спосіб заливання граней кольором, наприклад плавна зміна кольору від вершини до вершини, як це зроблено в наступному прикладі. Полігональний об'єкт створюється за допомогою функції patch. При його створенні використовуються властивості полігонального об'єкта:
· Vertices - матриця з трьох стовпців, у рядки якої записані координати вершин полігонального об'єкта;
· Faces - матриця, кожен рядок якої містить номера вершин, що утворять різної грані полігонального об'єкта
· FaceVertexCData - кольору вершин, у прикладі вони задаються у форматі RGB, колір установлюється зелений, а його відтінок залежить від висоти крапки над площиною xy.
· FaceColor - спосіб заливання граней полігонального об'єкта кольором, у прикладі використовується значення 'interp', тобто плавна зміна кольору від вершини до вершини.
При візуалізації опуклої оболонки множини крапок із тривимірного простору за допомогою функції patch необхідно самому встановлювати тривимірний вид осей, тому що patch відноситься до низькоурівневих граничних функцій, що не змінюють вид осей.
. Його початкові умови задовольняють умові 
(лежать на колі радіуса R у фазовій площині). Потрібно знайти значення 
, при яких енергія вихідного сигналу буде максимальної або мінімальною: 
.
Програму для керування процесом моделювання оформимо у вигляді окремого m_файлу optic.m.
(маятник відхиляємо вправо і штовхаємо від положення рівноваги). Мінімальна енергія вихідного сигналу,рівна Е = 1,264, досягається при j = 136про, тобто при початкових умовах 
маятник відхиляємо вліво і штовхаємо до положення рівноваги). Відповідні графіки приведені на рис. Б (вгоры).
). На рис. а (внизу) показані графіки енергії Е в залежності від кута _ і максимального і мінімального вихідного сигналу в у залежності від часу (рис., б) для R = 2,29, там же приведені чисельні значення екстремальних параметрів.
