Экономико-математическое моделирование : Диференційні рівняння як основа математичного опису енергетичної системи

Диференційні рівняння як основа математичного опису енергетичної системи

Контрольна робота з теми:

“Диференційні рівняння як основа математичного опису енергетичної системи”

1. Вихідні дані для реалізації системи звичайних диференційних рівнянь

Система диференційних рівнянь:

Початкові умови: А=0, В=1

t(0)=0, x(0)=0, y(0)=0

Задана точність: Е=

Обираємо с=6

1.1 Математична основа засобу Рунге - Кутта

Засіб Рунге -Кутта можливо получити, якщо разкласти у ряд Тейлора значення у(х)

y(x0+h)=y(x0)+h(x0)h3 +hnyn(x0)

xi=x(0)+Ih

yi+1=yi+•(K1i+2K2i+2K3i+2K4i)

K1i=h•f(xi,yi)

K2i=h•f(xi+•yi+)

K3i=h•f(xi+

K4i=h•f(xi+h•yi+K3)

Блок - схема головного модуля по Рунге - Кутту:

Реалізація програми за засобом Рунге - Кутта:

DECLARE SUB KUTT (T!, X!, Y!, A%, B%, C%, E!, H!, N%, T(), X(), Y(), K1X!, K1Y!, K2X!, K2Y!, K3X!, K3Y!, K4X!, K4Y!)

DECLARE SUB GRAF (T!, X!, Y!, A%, B%, C%, E!, H!, N%, T(), X(), Y(), K1X!, K1Y!, K2X!, K2Y!, K3X!, K3Y!, K4X!, K4Y!)

INPUT "C"; C%

E! = C% * 10 ^ (-4)

H! = E! ^ (1 / 4)

CONST A% = 0: CONST B% = 1

DIM SHARED T!(2000), X!(2000), Y!(2000), K1X!(2000), K1Y!(2000), K2X!(2000), K2Y!(2000), K3X!(2000), K3Y!(2000), K4X!(2000), K4Y!(2000)

T(0) = 0: X(0) = 0: Y(0) = 0

M1: CALL KUTT(T!, X!, Y!, A%, B%, C%, E!, H!, N%, T(), X(), Y(), K1X!, K1Y!, K2X!, K2Y!, K3X!, K3Y!, K4X!, K4Y!)

FOR I% = 0 TO N%

X1(I%) = X(I%)

Y1(I%) = Y(I%)

NEXT I%

H! = H! / 2

CALL KUTT(T!, X!, Y!, A%, B%, C%, E!, H!, N%, T(), X(), Y(), K1X!, K1Y!, K2X!, K2Y!, K3X!, K3Y!, K4X!, K4Y!)

FOR I% = 0 TO N%

IF ABS(X1(I%) - X(I%)) * (16 / 15) > E! THEN

GOTO M1

ELSE GOTO M2

END IF

IF ABS(Y1(I%) - Y(I%)) * (16 / 15) > E! THEN

GOTO M1

ELSE GOTO M2:

END IF

NEXT I%

M2: FOR I% = 1 TO N%

PRINT T(I%), X(I%), Y(I%)

NEXT I%

PRINT "H"; H!

INPUT K!

CALL GRAF(T!, X!, Y!, A%, B%, C%, E!, H!, N%, T(), X(), Y(), K1X!, K1Y!, K2X!, K2Y!, K3X!, K3Y!, K4X!, K4Y!)

END

SUB GRAF (T!, X!, Y!, A%, B%, C%, E!, H!, N%, T(), X(), X1(), K1X!, K1Y!, K2X!, K2Y!, K3X!, K3Y!, K4X!, K4Y!)

SCREEN 2

VIEW (170, 50)-(470, 150)

WINDOW (-1, 1.5)-(1, -1.5)

FOR I% = 0 TO N% - 1

PSET (T(I%), X(I%))

PSET (T(I%), Y(I%))

LINE (T(I%), X(I%))-(T(I% + 1), X(I% + 1))

LINE (T(I%), Y(I%))-(T(I% + 1), Y(I% + 1))

NEXT I%

LINE (-1, 0)-(1, 0)

LINE (0, -1.5)-(0, 1.5)

END SUB

SUB KUTT (T!, X!, Y!, A%, B%, C%, E!, H!, N%, T(), X(), Y(), K1X!, K1Y!, K2X!, K2Y!, K3X!, K3Y!, K4X!, K4Y!)

N% = (B% - A%) / H!

FOR I% = 0 TO N%

T(I%) = T(0) + I% * H!

K1X(I%) = H! * (-2 * X(I%) + 5 * Y(I%))

K1Y(I%) = H! * ((EXP(.5 * Y(I%) + T(I%)) - EXP(-.5 * Y(I%) + T(I%))) / 3 +.5 * Y(I%))

K2X(I%) = H! * (-2 * (X(I%) + K1X(I%) / 2) + 5 * (Y(I%) + K1Y(I%) / 2))

K2Y(I%) = H! * ((EXP(.5 * (Y(I%) + K1Y(I%) / 2) + (T(I%) + H! / 2) - EXP(-.5 * Y(I%) + K1Y(I%) / 2) - (T(I%) + H! / 2))) / 3 +.5 * (Y(I%) + K1Y(I%) / 2))

K3X(I%) = H! * (-2 * (X(I%) + K2X(I%) / 2) + 5 * (Y(I%) + K2Y(I%) / 2))

K3Y(I%) = H! * ((EXP(.5 * (Y(I%) + K2Y(I%) / 2) + (T(I%) + H! / 2) - EXP(-.5 * Y(I%) + K2Y(I%) / 2) - (T(I%) + H! / 2))) / 3 +.5 * (Y(I%) + K2Y(I%) / 2))

K4X(I%) = H! * (-2 * (X(I%) + K3X(I%)) + 5 * (Y(I%) + K3Y(I%)))

K4Y(I%) = H! * ((EXP(.5 * (Y(I%) + K3Y(I%)) + (T(I%) + H!) - EXP(-.5 * Y(I%) + K3Y(I%)) - (T(I%) + H!))) / 3 +.5 * (Y(I%) + K3Y(I%) / 2))

X(I% + 1) = X(I%) + 1 / 6 * (K1X(I%) + 2 * K2X(I%) + 2 * K3X(I%) + K4X(I%))

Y(I% + 1) = Y(I%) + 1 / 6 * (K1Y(I%) + 2 * K2Y(I%) + 2 * K3Y(I%) + K4Y(I%))

NEXT I%

END SUB

Результати реалізації системи диференційних рівнянь за засобом Рунге - Кутта

Нопт=0,07825423

1.2 Математична основа способу Мілна

Для реалізації засобу Мілна необхідно мати інформацію о попередніх точках. Тому засіб Мілна реалізуєтся після підрахунків по засобу Рунге-Кутта з заданной точностью.

Формула прогнозу:

Визначаємо значення проізводної:

Формула корекції:

Якщо, то закінчуємо розрахунок і даємо поманду на друк результату.

1.3 Реалізація програми за способом Мілна

DECLARE SUB MILN (T!, X!, Y!, A%, B%, H!, N%, E!, C%, X(), Y(), T(), LX3!, LY3!, LX2!, LY2!, LX1!, LY1!, XP!, YP!, XK!, YK!, MPX!, MPY!, MKX!, MKY!, XK1!, YK1!)

DECLARE SUB KUTT (T!, X!, Y!, A%, B%, H!, N%, E!, C%, X(), Y(), T(), KX1!, KY1!, KX2!, KY2!, KX3!, KY3!, KX4!, KY4!)

DECLARE SUB GRAF (T!, X!, Y!, A%, B%, H!, N%, E!, C%, X(), Y(), T(), LX3!, LY3!, LX2!, LY2!, LX1!, LY1!, XP!, YP!, XK!, YK!, MPX!, MPY!, MKX!, MKY!, XK1!, YK1!)

INPUT "C"; C%

E! = C% * 10 ^ (-4)

H! = E! ^ (1 / 4)

CONST A% = 0: CONST B% = 1

DIM SHARED T!(2000), X!(2000), Y!(2000), KX1!(2000), KY1!(2000), KX2!(2000), KY2!(2000), KX3!(2000), KY3!(2000), KX4!(2000), KY4!(2000)

DIM SHARED LX1!(2000), LY1!(2000), LX2!(2000), LY2!(2000), LX3!(2000), LY3!(2000), XP!(2000), YP!(2000), XK!(2000), YK!(2000), MPX!(2000), MPY!(2000), MKX!(2000), MKY!(2000), XK1!(2000), YK1!(2000)

T(0) = 0: X(0) = 0: Y(0) = 0

CALL KUTT(T!, X!, Y!, A%, B%, H!, N%, E!, C%, X(), Y(), T(), KX1!, KY1!, KX2!, KY2!, KX3!, KY3!, KX4!, KY4!)

FOR I% = 0 TO N%

PRINT T(I%), X(I%), Y(I)

NEXT I%

INPUT L!

CALL GRAF(T!, X!, Y!, A%, B%, H!, N%, E!, C%, X(), Y(), T(), LX3!, LY3!, LX2!, LY2!, LX1!, LY1!, XP!, YP!, XK!, YK!, MPX!, MPY!, MKX!, MKY!, XK1!, YK1!)

INPUT P!

CALL MILN(T!, X!, Y!, A%, B%, H!, N%, E!, C%, X(), Y(), T(), LX3!, LY3!, LX2!, LY2!, LX1!, LY1!, XP!, YP!, XK!, YK!, MPX!, MPY!, MKX!, MKY!, XK1!, YK1!)

FOR I% = 3 TO N%

IF (MPX(I% + 1) - MKX(I% + 1)) > E! THEN

XK1(I% + 1) = X(I% - 1) + (1 / 3) * H! * (MKX(I% + 1) + 4 * LX3(I%) + LX2(I% - 1))

PRINT T(I% + 1), XK1(I% + 1)

ELSE

XP(I% + 1) = XK(I% + 1)

PRINT T(I% + 1), XK(I% + 1)

END IF

FOR I% = 3 TO N%

IF (MPY(I% + 1) - MKY(I% + 1)) > E! THEN

YK1(I% + 1) = Y(I% - 1) + (1 / 3) * H! * (MKY(I% + 1) + 4 * LY3(I%) + LY2(I% - 1))

PRINT T(I% + 1), YK1(I% + 1)

ELSE

YP(I% + 1) = YK(I% + 1)

PRINT T(I% + 1), YK(I% + 1)

END IF

NEXT I%

INPUT M!

CALL GRAF(T!, X!, Y!, A%, B%, H!, N%, E!, C%, X(), Y(), T(), LX3!, LY3!, LX2!, LY2!, LX1!, LY1!, XP!, YP!, XK!, YK!, MPX!, MPY!, MKX!, MKY!, XK1!, YK1!)

END

SUB GRAF (T!, X!, Y!, A%, B%, H!, N%, E!, C%, X(), Y(), T(), LX3!, LY3!, LX2!, LY2!, LX1!, LY1!, XP!, YP!, XK!, YK!, MPX!, MPY!, MKX!, MKY!, XK1!, YK1!)

SCREEN 2

VIEW (170, 50)-(470, 150)

WINDOW (-1, 1)-(1, -1)

FOR I% = 0 TO N% - 1

PSET (T(I%), KX(I%))

PSET (T(I%), YK(I%))

LINE (T(I%), XK(I%))-(T(I% + 1), XK(I% + 1))

LINE (T(I%), YK(I%))-(T(I% + 1), YK(I% + 1))

NEXT I%

FOR I% = 4 TO N% - 1

PSET (T(I%), KX(I%))

PSET (T(I%), YK(I%))

LINE (T(I%), XK(I%))-(T(I% + 1), XK(I% + 1))

LINE (T(I%), YK(I%))-(T(I% + 1), YK(I% + 1))

PSET (T(I%), XK1(I%))

PSET (T(I%), YK1(I%))

LINE (T(I%), XK1(I%))-(T(I% + 1), XK1(I% + 1))

LINE (T(I%), YK1(I%))-(T(I% + 1), YK1(I% + 1))

NEXT I%

LINE (-1, 0)-(1, 0)

LINE (0, -1.5)-(0, 1.5)

END SUB

SUB KUTT (T!, X!, Y!, A%, B%, H!, N%, E!, C%, X(), Y(), T(), KX1!, KY1!, KX2!, KY2!, KX3!, KY3!, KX4!, KY4!)

N% = (B% - A%) / H!

I% = 0

DO

T(I%) = T(0) + I% * H!

KX1(I%) = H! * (-2 + X(I%) + 5 * Y(I%))

KY1(I%) = H! * ((EXP(.5 * Y(I%) + T(I%)) - EXP(-.5 * Y(I%) + T(I%))) / 3 +.5 * Y(I%))

KX2(I%) = H! * (-2 + (X(I%) + KX1(I%) / 2) + 5 * (Y(I%) + KX1(I%) / 2))

KY2(I%) = H! * ((EXP(.5 * (Y(I%) + KY1(I%) / 2) + (T(I%) + H! / 2) - EXP(-.5 * Y(I%) + KY1(I%) / 2 - (T(I%) + H! / 2))) / 3 +.5 * (Y(I%) + KY1(I%) / 2)))

KX3(I%) = H! * (-2 + (X(I%) + KX2(I%) / 2) + 5 * (Y(I%) + KY2(I%) / 2))

KY3(I%) = H! * ((EXP(.5 * (Y(I%) + KY2(I%) / 2) + (T(I%) + H! / 2) - EXP(-.5 * Y(I%) + KY2(I%) / 2 - (T(I%) + H! / 2))) / 3 +.5 * (Y(I%) + KY2(I%) / 2)))

KX4(I%) = H! * (-2 + (X(I%) + KX3(I%) / 2) + 5 * (Y(I%) + KY2(I%) / 2))

KY4(I%) = H! * ((EXP(.5 * (Y(I%) + KY3(I%) / 2) + (T(I%) + H! / 2) - EXP(-.5 * Y(I%) + KY3(I%) / 2 - (T(I%) + H! / 2))) / 3 +.5 * (Y(I%) + KY3(I%) / 2)))

X(I% + 1) = X(I%) + (1 / 6) * (KX1(I%) + 2 * KX2(I%) + 2 * KX3(I%) + KX4(I%))

Y(I% + 1) = Y(I%) + (1 / 6) * (KY1(I%) + 2 * KY2(I%) + 2 * KY3(I%) + KY4(I%))

I% = I% + 1

LOOP UNTIL I% > N%

END SUB

SUB MILN (T!, X!, Y!, A%, B%, H!, N%, E!, C%, X(), Y(), T(), LX3!, LY3!, LX2!, LY2!, LX1!, LY1!, XP!, YP!, XK!, YK!, MPX!, MPY!, MKX!, MKY!, XK1!, YK1!)

N% = (B% - A%) / H!

FOR I% = 3 TO N%

T(I%) = T(0) + I% * H!

LX3(I%) = -2 * X(I%) + 5 * Y(I%)

LY3(I) = (EXP(.5 * (Y(I%) + T(I%)) - EXP(-.5 * Y(I%) - T(I%))) / 3 +.5 * Y(I%))

LX2(I% - 1) = -2 * X(I% - 1) + 5 * Y(I% - 1)

LY2(I% - 1) = (EXP(.5 * (Y(I% - 1) + T(I% - 1)) - EXP(-.5 * Y(I% - 1) - T(I% - 1))) / 3 +.5 * Y(I% - 1))

LX1(I% - 2) = -2 * X(I% - 2) + 5 * Y(I% - 2)

LY1(I% - 2) = (EXP(.5 * (Y(I% - 2) + T(I% - 2)) - EXP(-.5 * Y(I% - 2) - T(I% - 2))) / 3 +.5 * Y(I% - 2))

XP(I% + 1) = X(I% - 3) + (4 / 3) * H! * (2 * LX3(I%) - LX2(I% - 1) + 2 * LX1(I% - 2))

YP(I% + 1) = Y(I% - 3) + (4 / 3) * H! * (2 * LY3(I%) - LY2(I% - 1) + 2 * LY1(I% - 2))

MPX(I% + 1) = 2 + XP(I% + 1) + 5 * YP(I% + 1)

MPY(I% + 1) = EXP(.5 * YP(I% + 1) + T(I% + 1)) - EXP(-.5 * YP(I% + 1) - T(I% + 1)) / 3 +.5 * YP(I% + 1)

XK(I% + 1) = X(I% - 1) + (1 / 3) * H! * (MPX(I% + 1) - 4 * LX3(I%) + LX2(I% - 1))

YK(I% + 1) = Y(I% - 1) + (1 / 3) * H! * (MPY(I% + 1) - LY3(I%) + LY2(I% - 1))

MKX(I% + 1) = -2 + XK(I% + 1) + 5 * YK(I% + 1)

MKY(I% + 1) = (EXP(.5 * YK(I% + 1) + T(I% + 1)) - EPX(-.5 * Y(I + 1) - T(I% + 1)) / 3 +.5 * YK(I% + 1))

NEXT I%

END SUB

1.4 Результати реалізації програми за способом Мілна

2. Реалізація контролю працездатності енергетичної системи

Початкові дані

Функціональний рівень зміни температури теплоносія, що гріє, для підігрівача гарячого водопостачання теплової підстанції (за варіантом №6)

Еталонний рівень зміни температури теплоносія, що гріє, 70….300С. Рівень підігріву місцевої води, що нагрівається, 5….600С.

1. Архітектура експертних систем

Розвиток та вдосконалювання обчислювальних комплексів, інформаційних технологій пов'язані з розробкою експертних систем, здатних обробляти не тільки кількісні дані, але й різного роду знання, проводячи аналіз поведінки енергетичних систем і приймаючи експертні рішення. Запропоновано керувати функціонуванням енергетичних систем на основі діагностичної інформації з використанням архітектури експертних систем, основою яких є динамічна система, що відбиває через характер реакцій на збурювання особливості функціонування енергетичних систем (її назва в експертній системі - динамічна підсистема). Іншими модулями, що входять до складу експертної системи, можуть бути блоки діагностування ситуації, ефективності, надійності тощо, з відповідним математичним описом і подальшим їх нарощуванням (рис. 1) [1-5].

Рис. 1 - Архітектура експертних систем: 1 - динамічна підсистема; 2 - модуль діагностування ефективності; 3 - модуль діагностування ситуації; 4 - модуль надійності (діагностування структурних параметрів)

2. Математичне моделювання енергетичної системи

Основою для контролю працездатності енергетичної системи є підігрівача гарячого водопостачання на тепловій підстанції. Ця основа здобута в результаті розв'язання системи нелінійних диференційних рівнянь передаточна функція по каналу "температура місцевої води, що нагрівається, - витрата теплоносія, що гріє" [2-4].

Температура роздільної стінки и:

Де С - питома теплоємність, Кдж/кг К;

D - витрата речовини, кг/с;

ТВ, ТМ - постійні часу, що характеризують теплову акумулюючу здатність робочого тіла, металу, с;

g - питома маса речовини, кг/м;

h - питома поверхня, м2/м;

t - температура робочого тіла, К;

z - координата довжини теплообмінника, м;

в - товщина стінки теплообмінника, м;

б - коефіцієнт тепловіддачі, кВт/м2К;

л - теплопровідність металу стінки теплообмінника, кВт/мК;

и, у - температура роздільної стінки, теплоносія, що гріє, К;

S (щj) - параметр перетворення Лапласа;

щ - частота.

Індекси:0 - стаціонарний режим;

1, 2 - вхід, вихід із теплообмінника;

в - потік робочого тіла;

н - потік теплоносія, що гріє;

м - металева стінка.

Для переходу із частотної області до області реального часу реалізую на ПЕОМ за способом Сімпсона такий інтеграл:

3. Логічне моделювання контролю працездатності

Контроль працездатності та ідентифікація стану енергетичної системи відбувається на основі графу причинно-наслідкових зв'язків динамічної підсистеми як основи експертної системи (рис. 2) [2-4].

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

Рис. 2 - Граф причинно-наслідкових зв'язків динамічної підсистеми: СТ - контроль події; Z - логічні відношення; ST - ідентифікація події. Індекси: 1 - Впливи; 2 - внутрішні параметри, що діагностуються; 3 - коефіцієнт рівнянь динаміки; 4 - істотні параметри, що діагностуються; 5 - динамічні параметри; С - контроль працездатності; S - стан

Повідомлення, що підтверджують ці властивості, одержувані динамічною системою від інших модулів експертної системи (якщо вони діагностуються), можуть бути використані для вироблення остаточного рішення та проведення подальших оперативних операцій [2-4].

Допустиму працездатність енергетичної системи визначаю такою логічною структурою:

. (6)

Рецепція ж такої результуючої інформації.

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

Індекси: рів. - рівень функціонування;

вст. - встановлення значення параметра;

розрах. - розрахункове значення параметра;

низ. - рівень зміни температури теплоносія, що гріє, 70…300С.