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



Сторінка3/5
Дата конвертації04.06.2017
Розмір0.69 Mb.
ТипМетодичні вказівки
1   2   3   4   5

X΄=F(t, X), X(t0)=X0. (3.6)

Змінні стану та їх початкові умови у процесі перетворення із форми (3.1) у форму (3.4) однозначно пов’язуються з невідомою функцією та її похідними.



Чисельний розв‘язок системи (3.4) полягає у приблизному визначенні для функцій , які є частковим розв‘язком системи ДР, таблиці її значень для деякої послідовності моментів часу на інтервалі . Якщо вектор часу сформувати у вигляді стовпця, то значення змінних стану утворять матрицю, як це показано у вигляді табл. 3.1. При цьому аналітичний розв‘язок залишається невідомим.

Основна ідея чисельних методів полягає у можливості розрахувати похідні від невідомих сигналів на кожному кроці чисельного інтегрування (ЧІ) і зробити цей крок у відомому напрямі від попередніх значень, у якості яких на першому кроці виступають початкові умови. У такий спосіб табл.2.1 заповнюється по рядкам зверху донизу. За розрахованими даними легко побудувати графіки перехідних процесів, як це показано на рис. 3.1.

Точність розв‘язку залежить від кроку ЧІ, тобто відстані між сусідніми інтервалами часу .

Серед чисельних методів розв‘язання ДР розрізняють одноступінчаті та багатоступінчаті.

В одноступінчатих методах для розрахунку значень змінних стану у точці , тобто. , необхідна інформація про одне попереднє значення розв‘язку . Одними з найпоширеніших одноступінчатих методів є методи Рунге-Кутти. Порядок методів Рунге-Кутта визначається кількістю розрахунків значень похідних на одному кроці ЧІ. Методи Рунге-Кутта можна використовувати як з постійним кроком, так і з автоматичним вибором кроку (АВК), який обирається з умов забезпечення за даної точності розв‘язку. Підвищити точність можна використанням комбінованих методів Рунге-Кутти різних порядків, які дають можливість визначити приблизно величину похибки і компенсувати її.

Багатоступінчаті методи для виконання одного кроку ЧІ потребують інформації про значення невідомих змінних у декількох попередніх точках часу, кількість яких визначає порядок методу. Серед багатокрокових методів виділяють велику групу під загальною назвою «методи прогнозу і корекції» (Predictor-Corrector Methods). Найбільш відомими серед них є методи Адамса, Бешфорта, Мілна, Гіра, Розенброка. Оскільки на першому кроці ЧІ відомими є тільки початкові умови, то методи прогнозу і корекції «запускаються» методами Рунге-Кутта відповідного порядку. При використанні методів прогнозу і корекції знайдений при першому наближенні розв‘язок потім уточнюється за допомогою ітераційних алгоритмів, в яких кількість ітерацій залежить від заданої точності розв‘язку.

При виборі певного методу розв‘язання ДР слід виходити з того, що одноступінчаті методи потребують менше оперативної пам‘яті, але працюють при заданій точності повільніше, ніж багатоступінчаті методи такого ж порядку. Також треба мати на увазі, що багатоступінчаті методи більш пристосовані для розв ‘язання так званих жорстких (stiff) ДР, для яких характерно сполучення швидкої та повільної динаміки, тобто відношення найбільшої сталої часу до найменшої з великим числом.




    1. Методи Рунге-Кутти

До числа найбільш розповсюджених чисельних методів рішення диференційних рівнянь належать методи Рунге-Кутти (Runge-Kutta methods).

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



(3.7)

Показник ступеню при в останньому підсумованому члені ряду Тейлора визначає порядок методу.

Метод Рунге-Куттипершого порядку називають методом Ейлера (Euler’s method), другого порядку – модифікованим методом Ейлера (improved Euler’s formula), або методом Ейлера-Коші, або Heun’s method, п‘ятого порядку – методом Дормана-Прінса (Dormand-Prince). Коли мова йде про метод Рунге-Кутта без визначення його порядку, то найчастіше мають на увазі метод Рунге-Кутти 4-го порядку, який застосовують найчастіше, завдяки досягненню компромісу між складністю, швидкодією та точністю методу.

Відповідно до метода Ейлера один крок рішення системи диференційних рівнянь (3.4) з початковими умовами (3.5) виконується за формулою:



(3.8)

Метод Ейлера-Коші потребує обчислення вектора похідних (правих частин диференційних рівнянь) F(t, Х) у двох точках:



(3.9)

(3.10)

Відповідно до використання метода Рунге-Кутта четвертого порядку вектор похідних на кожному кроці чисельного інтегрування обчислюється чотири рази :



(3.11)

. (3.12)

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

Похибка методів Рунге-Кутти визначається виразом:

.

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




    1. Розв‘язання систем лінійних диференційних рівнянь
      в пакеті MatLab


Серед методів розв‘язання ДР з фіксованим кроком в Simulink пропонуються методи Рунге-Кутта першого-п‘ятого порядків, а саме:

ode5метод Рунге-Кутта 5-го порядку (Dormand-Prince formula);

ode4метод Рунге-Кутта 4-го порядку (fourthe-order Runge-Kutta formula);

ode3метод Рунге-Кутта 3-го порядку (Bogacki-Shampine formula);

ode2модифікований метод Ейлера (improved Euler’s formula / Heun’s method);

ode1метод Ейлера (Euler’s method).

Серед методів розв‘язання ДР з заданою точністю, в Simulink і в MatLab пропонуються комбіновані методи Рунге-Кутти та деякі методи прогнозу та корекцій, а саме:



ode45комбінований метод Рунге-Кутта (4-5)-го порядку з автоматичним вибором кроку, призначений для розв‘язання нежорстких ДР;

ode23комбінований метод Рунге-Кутта (2-3)-го порядку з автоматичним вибором кроку, призначений для розв‘язання нежорстких та слабко жорстких ДР;

ode113багатокроковий метод Адамса-Моултона-Бешфорта (метод прогнозу та корекцій змінного порядку), призначений для розв‘язання нежорстких ДР; може бути більш ефективним, ніж ode45, при високій точності;

ode15s – багатокроковий метод (1-5)-го порядку (за замовчанням 5-го), що базується на формулах чисельного диференціювання (NDFs), але може використовувати і формули зворотного диференціювання (BDFs) у відповідності з методом Гіра, спрямований на розв‘язання жорстких ДР;

ode23s – однокроковий метод, що базується на модифікованій формулі Розенброка 2-го порядку; для деяких жорстких ДР виявляється при невисокій точності виявляється більш ефективним, ніж ode15s;

ode23t – метод трапецій з використанням „вільного” інтерполянта для розв‘язання помірно жорстких ДР;

ode23tb – комбінований метод, що використовує на першому етапі формулу трапецій, а на другому – зворотну диференційну формулу другого порядку; як і метод ode23s, може бути при невисокій точності більш ефективним, ніж ode15s, для розв‘язання жорстких ДР.

Функції ode1, …, ode5 можна обирати тільки при розрахунку перехідних процесів у середовищі Simulink, а усі інші із вище перелічених – як у середовищі Simulink, так і у програмах, написаних алгоритмічною мовою пакета MatLab.

Звернення до усіх цих функцій має однаковий формат (за виключенням імені функції). Розглянемо варіанти їх використання на прикладі функції ode23. У мінімальній конфігурації звернення до цієї функції виглядає так:

[t,X]=ode23(odeFun, Tspan, X0)

де

[t,X]вектори, які містять рішення (tвектор-стовбець часу, Xматриця вихідних змінних, в якій кількість рядків дорівнює кількості елементів у векторі часу, а кількість стовпчиків – порядку системи диференційних рівнянь n, тобто кількості невідомих змінних); i-й рядок матриці Х(i,:) містить у собі значення усіх невідомих змінних в i-й момент часу (i=1:length(t)), а j-й стовпчик матриці X (:,j) містить у собі значення однієї змінної у кожний момент часу (j=1:n);



Tspan=[t0 tf]інтервал часу, який визначається початковим t0 та кінцевим tf часами інтегрування (t0 – необов’язковий параметр, який за замовченням дорівнює 0, тоді Tspan=tf) або монотонно зростаючий чи монотонно спадаючий вектор часу: Tspan=[t0 t1 t2 … tf];

X0 – вектор початкових умов;

odeFun дескриптор m-функції (function handle), у якій міститься опис системи диференційних рівнянь, тобто формули для розрахунку вектору перших похідних від невідомих змінних, або, як їх називають, функції правих частин системи диференційних рівнянь у нормальні й формі Коші (див. рівняння (10.5)) у заданий момент часу t.

Дескриптор функції конструюється символом @, за яким прямує ім‘я функції: @FunName. Замість нього можна використати ім‘я функції у вигляді строки символів ‘FunName’.

Обов‘язкова частина функції правих частин має такий формат:

function dX=MainRight (t,X)

dX(1)=…;


dX(2)=…;

dX(n)=…;



dX=dX(:);

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

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

» odeset


AbsTol: [ positive scalar or vector {1e-6} ]

BDF: [ on | {off} ]

Events: [ on | {off} ]

InitialStep: [ positive scalar ]

Jacobian: [ on | {off} ]

JConstant: [ on | {off} ]

JPattern: [ on | {off} ]

Mass: [ {none} | M | M(t) | M(t,y) ]

MassSingular: [ yes | no | {maybe} ]

MaxOrder: [ 1 | 2 | 3 | 4 | {5} ]

MaxStep: [ positive scalar ]

NormControl: [ on | {off} ]

OutputFcn: [ string ]

OutputSel: [ vector of integers ]

Refine: [ positive integer ]

RelTol: [ positive scalar {1e-3} ]

Stats: [ on | {off} ]

Vectorized: [ on | {off} ]

Значення опцій за замовчанням подані у фігурних дужках. Найбільший інтерес мають параметри AbsTol (абсолютна похибка), InitialStep (початковий крок), MaxStep (максимальний крок), RelTol (відносна похибка).

Змінити значення опцій можна за допомогою присвоєння вихідному параметру options значення функції odeset з чергуванням вхідних параметрів «ім‘я властивості» та «значення властивості» і доповнення виклику функції ode... сформованим у такий спосіб вхідним параметром. Наприклад, наведена нижче послідовність операторів

opt = odeset('RelTol',1e-4,'AbsTol',[1e-5 1e-5 1e-5], ‘MaxStep’, 0.01);

[t,y]=ode113(@mydiffur, 1.2, zeros (1,3), opt);

розв‘язує систему ДР 3-го порядку, описану у функції mydiffur, з нульовими початковими умовами на інтервалі часу методом Адамса-Моултона-Бешфорта з відносною точністю 10-4, абсолютною точністю по кожній змінній стану 10-5 і максимальним кроком 0.01.

Інформацію про поточні значення опцій після їх зміни можна побачити за допомогою оператора odeget, наприклад:

» odeget(opt,'RelTol')

ans =


0.0001

    1. Приклад розрахунку перехідних процесів в електричному лінійному колі в пакеті MatLab

      1. Умова задачі

Математичний опис лінійних електричних кіл створюють, використовуючи закони Ома та Кірхгофа з урахуванням тієї обставини, що напруги і струми на резисторі індуктивності та конденсаторі пов‘язані між собою співвідношеннями:

, , . (3.13)

Як приклад, розглянемо розгалужену схему, зображену на рис. 2.7, для якої треба розрахувати перехідні процеси, тобто знайти i1(t),
Рисунок 3.2

i2(t), i3(t), UC(t), при замиканні ключа при таких параметрах: E = 100 B,
J
k = 5 A, C = 200 мкФ, L = 0,3 Гн, r1 = 20 Om, r2 = 50 Om, r3 = 70 Om.

      1. Математичний опис

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

, (3.14)

, (3.15)

, (3.16)

. (3.17)

Рівняння (3.15), (3.16) – диференційні; (3.14), (3.17 ) – алгебраїчні рівняння, які часто називають рівняннями зв’язку.

Для використання чисельних методів розв‘язання ДР перетворимо їх у нормальну форму Коші:

З (3.18)

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

Об‘єднаємо невідомі змінні, від яких визначено похідні, в один вектор:

(3.19)

Ці сигнали утворюють вектор змінних стану. Їх кількість визначає порядок системи ДР.

З рівнянь зв‘язку визначаємо інші невідомі сигнали через змінні стану:

(3.20)

Розрахуємо початкові умови для заданої схеми, які збігаються з усталеними значеннями відповідних сигналів при розімкненому стані ключа. Така схема зображена на рис. 3.3.



Рисунок 3.3

В усталеному режимі усі похідні дорівнюють 0, звідси

(3.21)


      1. Розрахунок перехідних процесів програмним методом

Для розрахунку перехідних процесів програмним методом за допомогою m-функцій, розглянутих у попередньому підрозділі, розробимо на основі виразів (3.18)-(3.20) функцію, що буде обчислювати праві частини ДР (3.18), тобто вектор похідних від змінних стану. При цьому треба дотримуватися таких правил:

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

  • присвоїмо невідомим змінним у скалярному вигляді значення індексованих змінних стану (див. (3.19));

  • напишемо рівняння зв‘язку (3.20) у скалярній формі;

  • напишемо диференційні рівняння (3.19);

  • перетворимо вектор похідних у стовбець.

Після кожного з операторів присвоєння ставимо «крапку із зап‘ятою, щоб запобігти виведення зайвої проміжної інформації на екран.

% Математичний опис електричного кола

function dx=right1(t,x)

global E Jk R L C % Глобальні зміні

% Вираз невідомих змінних через змінні стану в скалярній формі

Uc=x(1); i3=x(2);

% Визначення допоміжних змінних із рівнянь зв‘язку в скал. формі

i1=(E+y(1))/R(1); i2=y(2)-Jk;

% Диференційні рівняння в нормальній формі Коші

dx(1)=(i2-i1)/C;

dx(2)=(E-i1*R(1)-i2*R(2)-y(2)*R(3))/L;

dx=dx(:); %Витягування вектору похідних у стовбець

Створену функцію слід зберегти у файлі, ім‘я якого співпадає з ім‘ям функції, тобто right1.m.

Після цього розробляємо головну програму (script file). При цьому будемо дотримуватися такої послідовності:



  • імена змінних, значення яких потребують як головна програма, так і функція правих частин, описуємо як глобальні, після чого визначаємо вихідні данні та розраховуємо початкові умови;

  • задаємось інтервалом часу розрахунку, який може бути уточнений після декількох спроб;

  • розв‘язуємо ДР, визначені у функції одним із методів, описаних у попередньому підрозділі; у такий спосіб визначаємо вектор-стовбець часу та матрицю змінних стану;

  • присвоїмо векторам невідомих змінних ДР значення відповідних стовбців матриці змінних стану;

  • розрахуємо з рівнянь зв‘язку (2.39) вектори допоміжних змінних;

  • на підставі отриманих чисельних даних побудуємо графіки перехідних процесів.

Після операторів, що визначають вихідні дані, «крапку із зап‘ятою» можна не ставити, а після операторів, що розраховують перехідні процеси, «крапку із зап‘ятою» краще поставити, щоб запобігти виведенню на екран значних масивів даних, які використовуються для побудови графіків.

% Текст головної програми

global E Jk R L C % Глобальні змінні

E=100, Jk=5, C=200e-6, L=0.3, R=[20 50 70] % Вихідні данні

X0=[-Jk*R(3) Jk] % Початкові умови

tk=0.1 % Час перебігу перехідних процесів

[t,X]=ode23 (@right1, tk, X0); % Розв‘язання системи ДР

Uc=X(:,1); I3=X(:,2); % Невідомі змінні ДР

I1=(X(:,1)+E)/R(1); I2=X(:,2)-Jk; Ic=I2-I1; % Допоміжні змінні

% Побудова графиків перехідних процесів

figure, subplot(121), plot(t,[I1 I2 I3 Ic]), grid on, legend('i1','i2','i3','ic');

subplot(122), plot(t,Uc), grid on

Результат роботи програми наведено на рис. 3.3.

Рисунок 3.3

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

Якщо ЕРС та/або струми джерел змінюються у часі, то закон зміни цих сигналів треба описати у функції правих частин у скалярному вигляді, а в головній програмі – у векторному.

Наприклад, якщо

, , ,

то у головній програмі зміняться перші 2 рядки

global Em w Jk R L C % Глобальні змінні

Em=220, f=50, w=2*pi*f, Jk=5, C=200e-6, L=0.3, R=[20 50 70]

та перед розрахунком допоміжних змінних приєднається оператор

Е=Em*sin (w*t);

Точно такий же оператор додаємо і у функцію правих частин ДР.

Для формування періодичних сигналів прямокутної та трикутної форми, у тому числі і пилкоподібної, у Signal Processing Toolbox існують функції square та sawtooth.

Функція square генерує прямокутний періодичний сигнал одиничної амплітуди з періодом 2*pi. Необов‘язковий параметр duty cycle визначає у процентах долю періоду, коли вихідний сигнал має додатне значення. За замовченням значення цього параметру становить 50%.

Функція sawtooth генерує трикутний періодичний сигнал одиничної амплітуди з періодом 2*pi. Необов‘язковий параметр width визначає долю періоду, коли вихідний сигнал досягає максимального значення. За замовченням значення цього параметру становить 1. При цьому вихідний сигнал має пилкоподібну форму з крутим заднім фронтом. При width =0 він стає пилкоподібним із крутим переднім фронтом, а при width =0.5 – набуває форму трикутної симетричної хвилі.

Для генерації перших 10 періодів пилкоподібного та прямокутного періодичних сигналів з частотою 50 Hz із кроком 1/10 kHz, можна скласти таку програму:

fs = 10000;

t = 0:1/fs:0.2;

x1 = sawtooth(2*pi*50*t);

x2 = square(2*pi*50*t);

subplot(211),plot(t,x1), axis([0 0.2 -1.2 1.2])

xlabel('Time (sec)');ylabel('Amplitude'); title('Sawtooth Periodic Wave')

subplot(212),plot(t,x2), axis([0 0.2 -1.2 1.2])

xlabel('Time (sec)');ylabel('Amplitude'); title('Square Periodic Wave')

Результати її виконання відображені на рис. 3.4.

Рисунок 3.4



Для формування одиночних (неперіодичних) трикутних, прямокутних та гаусовських імпульсів у Signal Processing Toolbox існують функції tripuls та rectpuls відповідно.

Імпульси формуються симетричними відносно точки . Ширина імпульсів за замовчанням дорівнює 1. Її можна змінити додатковим параметром.

Зразок використання цих функцій наведено нижче:

fs = 10000;

t = -0.1:1/fs:0.1;

x1 = tripuls (t, 20e-3);

x2 = rectpuls (t, 20e-3);

subplot (211), plot (t, x1), ylim ([-0.2 1.2])

xlabel('Time (sec)'); ylabel('Amplitude'); title('Triangular Aperiodic Pulse')

subplot (212), plot (t, x2), ylim ([-0.2 1.2])

xlabel ('Time (sec)'); ylabel ('Amplitude'); title ('Rect Aperiodic Pulse')

Їх вигляд показано на рис. 3.5



Рисунок 3.5





      1. Поділіться з Вашими друзьями:
1   2   3   4   5


База даних захищена авторським правом ©divovo.in.ua 2017
звернутися до адміністрації

войти | регистрация
    Головна сторінка


загрузить материал