Міністерство освіти і науки України
Національний університет “Львівська політехніка”
Лабораторна робота № 3
з курсу “Моделювання процесів
та елементів систем керування”
«Прискорений пошук усталених режимів
електромаґнетних елементів систем керування»
Мета роботи – вивчити методи прискореного пошуку усталених режимів електромаґнетних елементів систем керування, а також дослідити ці режими для заданої схеми використовуючи числові методи розв’язування систем нелінійних диференціальних рівнянь.
ОСНОВНІ ТЕОРЕТИЧНІ ВІДОМОСТІ
Таким чином, маючи систему диференціальних рівнянь, що описує стан пристрою і початкові умови, з допомогою чисельних методів ми можемо розрахувати перехідні процеси що в них протікають. Наступною задачею, яка викликає великий інтерес є аналіз усталених режимів роботи електромаґнетних елементів систем керування. Рівняння динаміки дають можливість проводити аналіз перехідних процесів. Інтеґруючи ці рівняння на достатньо великому інтервалі часу, отримаємо усталений режим. Безумовно, що такий підхід має два недоліки: 1) накопичення похибок чисельного інтеґрування; 2) значні затрати машинного часу. Тому розглянемо ряд методів, які позбавлені цих недоліків, дають можливість отримувати розв’язок в часовій області і з наперед заданою точністю збіжності ітераційного процесу. Як і усі ітераційні методи вони дуже зручні при програмуванні.
1.1. Модель чутливостей до початкових умов
Один з методів, що дає можливість проводити аналіз усталених режимів базується на ітераційних циклах Ньютона, або як його ще називають – модель чутливостей до початкових умов. Він дає можливість знаходити такі початкові умови, які при інтеґруванні рівнянь динаміки приводять до усталеного режиму.
В усталеному режимі змінні стану повинні задовольняти умову періодичності, яку запишемо у вигляді певного нелінійного рівняння
(1)
де – період вхідної дії.
Нелінійне рівняння (1) розв'язуємо ітераційним методом Ньютона
(2)
де – якобіан рівняння періодичності. Диференціюючи (2) за отримаємо
(3)
де – матриця монодромії. Запишемо її у вигляді добутку двох інших матриць. Для цього часткову похідну (3) запишемо як похідну складної функції
(4)
де – матриця коефіцієнтів рівнянь динаміки; – матриця додаткових чутливостей. Останню визначаємо в результаті інтеґрування додаткової системи диференціальних рівнянь першої варіації. Для цього продиференціюємо за часом матрицю S
. (5)
Запишемо диференціальні рівняння прикладу розглянутого в лабораторній роботі № 2
, (6)
де – вектор змінних стану; – матриця коефіцієнтів
. (7)
Тут – матриці коефіцієнтів, причому
Запишемо рівняння (5) як похідну складної функції
. (8)
Обчислимо часткову похідну для матриці коефіцієнтів
(9)
Враховуючи, що , , матрицю (9) запишемо у вигляді
(10)
Маючи рівняння динаміки можна обчислити часткові похідні матриці , а саме
, (11)
де , .
Підставивши (11) в (10) отримаємо остаточний вираз для матриці
(12)
Рівняння (5) запишемо у матричній формі
(13)
Рівняння (13) завжди лінійне. Його інтеґруємо сумісно з нелінійним рівнянням (6). Початкове наближення задаємо довільним, як правило, нульовим. Початкові умови згідно (2.6), (2.7) повинні строго задовольняти умову
(14)
Ітераційний процес (2) зупиняємо при виконанні нерівності
(15)
де – помилка збіжності ітераційного процесу.
1.2. Екстраполяційний метод
Як було показано вище, прискорений пошук вимушених періодичних режимів зводиться до обчислення початкових умов, які дають можливість в процесі інтеґрування рівнянь стану системи отримати безпосередньо усталений режим в обхід перехідного. Застосування методу Ньютона (2) для цієї цілі передбачає обчислення матриці монодромії.
Визначення початкових умов, що виключають перехідну реакцію системи, можна проводити екстраполяційним методом, який не залежить від природи диференціальних рівнянь кола і тому позбавлений процедури обчислення матриці монодромії.
Екстраполяційні методи цільову функцію використовують у вигляді Інтеґруючи (2.3) на періодах, породжуємо послідовність значень
(16)
де – період вхідної дії. Для послідовності (16), починаючи з , застосовуємо екстраполяцйну формулу
, (17)
де – початкові умови входження в усталений режим.
В якості функції EXTR доцільно використати – алґоритм, який виконує обчислення границі послідовності з експоненційними складовими. Формула для обчислення наступного значення має вигляд
(18)
де Результат екстраполяції, що відповідає EXTR в (17) рівний В (18) використовується процедура обертання Самельсона де – k-й елемент n-мірної колонки .
Для систем розмірності значення . У нашому прикладі =3. Початкові умови усталеного режиму для швидкозатухаючих перехідних реакцій компонент визначаємо інтеґруванням рівнянь стану системи на періодах. Як правило . На жаль, не існує строгого критерію вибору i , тому тут можливий лише евристичний підхід. Основний недолік екстраполяційних методів полягає в необхідності інтеґрування рівнянь динаміки на значному інтервалі часу.
Алґоритм обчислень
1. Інтеґруємо рівняння (6), від заданих початкових умов на періодах і визначаємо початкові умови періодичного режиму швидкозатухаючих компонент .
2. Маючи на -й ітерації початкові умови змінних стану (на першій ітерації умови п.1), інтеґруємо рівняння (6) на 6 періодах і породжуємо послідовність
. (19)
3. Згідно (18) визначаємо уточнене значення початкових умов
(20)
4. Перевіряємо умову (15) збіжності ітераційного процесу. Якщо вона не виконується, то процес повторюємо з п.2, в противному випадку зупиняємо ітераційний процес.
Текст основної програми
#include "stdafx.h"
#include "math.h"
#include "fstream"
#include "iostream"
#include "conio.h"
using namespace std;
const double a1=190,a2=230,a3=230,m1=3,m2=210,r1=5,r2=7,r3=11,C1=0.008,C2=0.0007,C3=0.003,Um=330,Rn2=19,Rn3=12,i1=0.45,i2=35,R2=10;
const double h=0.0001;
double ff (double psi)
{
if (psi<=0.4)
return 2*psi;
if ((psi<1.2)&&(psi>0.4))
return (1.2+(-2.375)*psi+(-0.625)*psi*psi+10.15625*psi*psi*psi);
if (psi>1.2)
return 40*psi+15-40*1.2;
}
double aii (double psi)
{
if (psi<=0.4)
return 2;
if ((psi<1.2)&&(psi>0.4))
return (-2.375+(-1.25)*psi+30.46875*psi*psi);
if (psi>1.2)
return 40;
}
void DfDt (double K[5], double X[5], double t)
{
double g1,g2,g3,g;
g=a1+a2+a3+aii(X[0]);
g1=a1/g;
g2=a2/g;
g3=-a3/g;
double B[5][5]={{g1,g2,g3,0,0},{-a2*g1, a2*(1-g2),-a2*g3,0,0},{a3*g1,a3*g2,a3*(1+g2),0,0},{0,0,0,1/C2,0},{0,0,0,0,1/C3}};
double I[3]={ff(X[0]-X[1]+X[2]),X[1],X[2]}, U[3]={Um*sin(2*3.14159*50*t),-X[3],-X[4]},
R[3][3]={{r1,0,0},{0,r2+Rn2,R2},{0,R2,r3+Rn3}}, Z[3]={0,0,0};
for (int i=0;i<3;i++)
for (int j=0;j<3;j++)
Z[i]+=R[i][j]*I[j];
for (int i=0;i<3;i++)
Z[i]=U[i]-Z[i];
double Y[5]={Z[0],Z[1],Z[2],X[1],X[2]};
for (int i=0;i<5;i++)
{
K[i]=0;
for (int j=0;j<5;j++)
K[i]+=B[i][j]*Y[j];
}
}
void rungecutta (double X[5])
{
fstream file;
file.open("lab3.dat", ios_base::out|ios_base::trunc);
double K1[5],K2[5],K3[5],K4[5],Z[5];
for (double t=0; t<1.0/50;t+=h)
{
DfDt(K1,X,t);
for (int i=0;i<5;i++) Z[i]=X[i]+h/2*K1[i];
DfDt(K2,Z,t+h/2);
for (int i=0;i<5;i++) Z[i]=X[i]+h/2*K2[i];
DfDt(K3,Z,t+h/2);
for (int i=0;i<5;i++) Z[i]=X[i]+h*K3[i];
DfDt(K4,Z,t+h);
for (int i=0;i<5;i++) X[i]=X[i]+h/6*(K1[i]+2*K2[i]+2*K3[i]+K4[i]);
cout<<t<<" "<<X[0]<<" "<<X[1]<<" "<<X[2]<<" "<<X[3]<<" "<<X[4]<<endl;
file<<t<<" "<<X[0]<<" "<<X[1]<<" "<<X[2]<<" "<<X[3]<<" "<<X[4]<<endl;
}
file.close();
}
void _tmain(int argc, _TCHAR* argv[])
{
const int m=5, q=5, n=11, p=2;
double X[5]={0,0,0,0,0},e[12][12][5], Eps=0.000001;
double V[5], sum;
bool cond;
int count=0;
for (int j=0;j<=n;j++)
for (int i=0;i<m;i++)
e[j][0][i]=0;
do
{
cond=false;
for (int j=1;j<=p;j++)
rungecutta (X);
for (int i=0;i<m;i++)
e[0][1][i]=X[i];
for (int j=0;j<10;j++)
{
rungecutta (X);
for (int i=0;i<m;i++)
e[j+1][1][i]=X[i];
if (j==0)
{
for (int i=0;i<=m-1;i++)
cond=cond||(fabs((e[1][1][i]-e[0][1][i])/e[1][1][i])*100>Eps);
if (cond==false) break;
}
}
if (cond)
{
for (int k=1;k<=n-1;k++)
for (int j=0;j<=n-k-1;j++)
{
for (int i=0;i<=m-1;i++)
V[i]=e[j+1][k][i]-e[j][k][i];
sum=0;
for (int i=0;i<=m-1;i++)
sum=sum+V[i]*V[i];
for (int i=0;i<=m-1;i++)
V[i]=V[i]/sum;
for (int i=0;i<=m-1;i++)
e[j][k+1][i]=e[j+1][k-1][i]+V[i];
}
for (int i=0;i<=m-1;i++)
X[i]=e[0][n][i];
count++;
}
}
while (cond);
cout<<count;
getch();
}
Результати виконання роботи
Висновок – отже, ми вивчили методи прискореного пошуку усталених режимів електромаґнетних елементів систем керування, а також дослідили ці режими для заданої схеми використовуючи числові методи розв’язування систем нелінійних диференціальних рівнянь.