CC BY
3Электронный журнал «Труды МАИ». Выпуск № 50
www.mai.ru/science/trudy/
УДК 621.372.061
Моделирование распределенных и дробно-распределенных процессов и систем управления в СКМ спектральным методом
В.В. Рыбин
Аннотация
В настоящее время для моделирования систем управления спектральным методом применяются различные системы компьютерной математики (СКМ) и пакеты их расширения. Современные технологии, связанные с фрактальным подходом в различных прикладных областях, в частности в теории динамических систем, порождают новую элементную базу, математические модели которой содержат дробные интегрирующие и дифференцирующие операторы.
В данной статье спектральный метод применяется для моделирования распределенных и дробно-распределенных процессов и систем управления. Для моделирования таких систем модифицирован пакет расширения МЬ8У_БМ СКМ МаШсаё. Применение этого пакета демонстрируется на примере решения уравнения аномальной диффузии и моделировании системы управления ядерной энергетической установкой с распределенными параметрами, которая содержит дробный ПИД регулятор.
Ключевые слова
нестационарные системы автоматического управления; спектральная форма математического описания; системы компьютерной математики; дробные интегрирующие и дифференцирующие звенья; аномальная диффузия; ядерная энергетическая установка.
1. Введение
В настоящее время теория дробных операторов находит все большее применение в теории управления и других предметных областях [1-4]. Микро- и нанотехнологии позволяют создавать технические элементы, которые физически реализуют дробные интегральные и дифференциальные операторы. Для систем автоматического управления предложена методика проектирования ПИД регуляторов дробного порядка [4]. В работе [5] спектральный метод [6-12] развить на нестационарные системы управления, содержащие дробные интегрирующие и дифференцирующие звенья, а для моделирования дробных систем управления летательными аппаратами модернизированы пакеты
расширения MLSY_SM, Spektr_SM+Simulink+Matlab, Spektr_SM+VisSim+Mathcad СКМ [12-18] расчета нестационарных непрерывных, дискретных и непрерывно-дискретных систем управления летательными аппаратами с сосредоточенными параметрами.
Заметим также, что для решения современных энергоемких задач работы космических аппаратов требуется энергия, дать которую в настоящее время способны только ядерные энергетические установки (ЯЭУ). Модели ЯЭУ описываются диффузионными уравнениями, а системы управления могут содержать дробные ПИД регуляторы. Спектральный метод удобен для моделирования таких систем управления ЯЭУ. Поэтому для моделирования таких систем управления модифицирован пакет расширения MLSY_SM СКМ Mathcad.
В данной статье спектральный метод применяется для моделирования распределенных и дробно-распределенных процессов и систем управления с использованием пакета расширения MLSY_SM СКМ Mathcad. Применение этого пакета демонстрируется на примере решения уравнения аномальной диффузии [20] и моделировании системы управления ядерной энергетической установкой с распределенными параметрами [21], которая содержит дробный ПИД регулятор.
2. Основные характеристики спектральной формы математического описания одномерных распределенных систем управления
Спектральные характеристики для представления функций времени и/или вектора состояния, а также линейных операторов определяются с использованием табличной формы представления многомерных матриц [9,10]. Рассмотрим здесь только некоторые спектральные характеристики описания одномерных распределенных процессов и систем, определенных относительно ортонормированных систем функций.
Пусть (в, x) е Q, где в еТ = [ü, t] - конечный промежуток времени,
x е Q = [a,b] ^ R, а система функций \e(i0,\,в,x)}^ 0 образует ортонормированный базис
пространства L2 (Q). Тогда гиперстолбец H(2,0) = (h^ ), элементы которой представляют
собой коэффициенты разложения функции h(e, x) е L2 (Q) в ряд по функциям базисной системы, называется спектральной характеристикой функции времени и координаты состояния h(e,x), т.е.
S [h(e, x)] = H (2,0) (1)
где
hVl =(e(io,ii,e,x),h(e,x))Li(Q), io,i, = 0,1,2,.... (2)
Обратный переход от спектральной характеристики к соответствующей функции времени осуществляется по формуле обращения
к(в,х) = ^1[Я(2,0)] = ]Г]Гк^ •в(10,'^в,х), (в,х) е 0 . (3)
'0
'„ =0 и =0
Заметим, что если система функций {з(',в)}^0 образует базис пространства (Г), а система {р^, х)}"=0 является базисом пространства (О) и к(в, х) = к0 (в) • к (х), где к0 (в) е (Г), к (х) е (О), то спектральную характеристику функции к(в, х) , определенную относительно базисной системы {е('0,' ,в, х) = ^('0 ,в) • р^, х)}^ ,
можно представить в виде
Н (2,0) = Н0 (1,0) ® Н (1,0), (4)
где Н (1,0) - спектральная характеристика функции к0 (в), определенная относительно базисной системы {р(г0 ,в)}^=0; Н (1,0) - спектральная характеристика функции к (х),
определенная относительно базисной системы {р(\, х)}"=0.
В дальнейшем, по мере необходимости, при описании систем управления в спектральной области, спектральные характеристики (1) будем называть гиперстолбцовыми матрицами нестационарных спектральных характеристик (ГСМ НСХ).
Линейные операторы. Пусть М - линейный оператор, определенный на пространстве (0), а {е('0,',в,х)}^ 0 - ортонормированный базис этого пространства,
б = Т хО, веТ = [0, х е О = [а, Ь] с Я.
Гиперквадратная матрица М(2,2) = (тго^оЛ), элементы которой определяются
формулой тШоЛ = (e(io,/1,в,x),Me(jo,jl,в,х))^(ау io,КЛ^К = ^А---> называется спектральной характеристикой оператора М, т.е.
з|м]=М(2,2) . (5)
Заметим, что
эМкв, х)]=М (2,2) • Н (2,0), (6)
т.е. спектральная характеристика образа функции к(/, х) равна произведению
спектральной характеристики оператора М и спектральной характеристики функции к(Г, х).
Приведем свойства спектрального преобразования линейных операторов, определенных на пространстве Ь2 (0).
1. Спектральное преобразование тождественного оператора. Спектральная характеристика тождественного оператора I представляет собой единичную матрицу размерности 4 : Б[1] = Е(2,2) .
2. Спектральное преобразование композиции операторов. Спектральная характеристика композиции линейных операторов М и N равна произведению их спектральных характеристик:
б[М ° N ] = Б[М ]• ]=М (2, 2) • N (2,2), (7)
где М (2, 2) = Б|ММ ], N (2, 2) = ].
3. Спектральное преобразование обратного оператора. Предположим, что для
линейного оператора М существует обратный оператор М-1. Тогда спектральная характеристика обратного оператора равна обратной спектральной характеристике оператора М~ :
б|М 1 ]=М-\2,2), (8)
где М (2,2) = б|м ].
Примеры линейных операторов
1. Операторы умножения. Линейный оператор А называется оператором умножения на функцию а(в, х), если для любой функции И(в, х) е (@) справедливо
выражение АН(в, х) = а(в, х)Н(в, х), где а(в, х) - локально интегрируемая функция на множестве Q = Т хй, в еТ = [0, ^ ], х ей = [а, Ь] ^ Я. Тогда гиперквадратная матрица А(2,2)
= (о^КоЛ), элементы которой определяются формулой
а н =(е(*о,Ь,в,х),а(в,х)• е(г0,^,в,х))^в), /0,,К,К = 0,1,2,... называется спектральной
характеристикой оператора умножения на функцию а(в, х), т.е. б[а] = А(2,2) .
Заметим, что если система функций {<р(г,в)}°°=0 образует базис пространства £2 (Т), а система {р(\, х)}"=0 является базисом пространства £2 (й) и а(в, х) = а0 (в) • а (х), где а0 (в) е £2 (Т), а (х) е £2 (й), то спектральную характеристику оператора умножения а(в,х), определенную относительно базисной системы {е(/0,\,в,х) = <р(г0,в) • р(\,х) . = 0,
можно представить в виде
А(2,2) = Ао(1,1) 0 4(1,1) (9)
где А0 (1,1) - спектральная характеристика оператора умножения на функцию а0(в),
определенная относительно базисной системы {<р(г,в)}°°=0; А (1,1) - спектральная
характеристика оператора умножения на функции а (х) , определенная относительно
базисной системы {р(\, х)}"=0.
2. Операторы интегрирования. Рассмотрим операторы интегрирования, заданные на пространстве ¿2 (0). Будем рассматривать оператор дробного интегрирования по времени [1, 5]:
Щк(0, х) = Г(в~ТУ'1 к(т, х)йт, (10)
1 0 Г(/)
а также оператор дробного интегрирования по координате состояния [1, 5]:
х) = Г( к(в, /Щ. (11)
1 а г/)
Эти операторы при / = 1,2,...,т превращаются в обычные интегралы т -го порядка [1].
Спектральные характеристики операторов дробного интегрирования и
, определенные относительно базисной системы {е('0,',в,х)}°°.=0 пространства Х2 (б), вычисляются по определению спектральных характеристик линейных операторов и обозначаются (2,2) = (Р " ), (2,2) = (Р_" ) соответственно:
4 У V '0170 Л ' 1 1г1'2 Л72У
=
'0.'1.>0.Л
'0 «0 71 ' 1 1г1'2
^ в (0 — тТ_! ^
о Г(/) ^(е)
Р^/70 71 = (^('0, '1,в, х), Г ^ /У1 1 р(/0, о^л а г(/)
, i0,'1^.70^.71 = 0,l,2,..
^ ¿2(2)
В случае если функции базисной системы пространства ¿2 (б) представляются в
виде {е('0,',в,х) = ^('0,в)• р(',х)}° ,=0, т.е. в виде произведения функций базисных систем
пространств ¿2 (Т) и ¿2 (О), то спектральные характеристики операторов дробного интегрирования представляются в виде
Р(2,2) = Р_"(1,1) ® £(1,1), Р;м (2,2) = £(1,1) ® Р л (1,1). (12)
Здесь двумерная матрица Р "(1,1) = (Р/"), элементы которой определяются формулой
( в.(а_ т\"_1 Л
Р. _" =
У
<р(',в), Г{в_У" т)ёт , ',7 = 0,1,2,., (13)
Г Г(") Л2(т)
называется спектральной характеристикой оператора интегрирования дробного порядка М> 0 по времени, т.е. б[оо+^]= P~м(1,1) [5], а двумерная матрица Р~м(1,1) = (Рьм), элементы которой определяются формулой
( р^-1 Л
р-м — Р1К
p(i, х), |(х-РМ , Ч, к = 0,1,2,..., (14)
(м) )ь2 (й)
называется спектральной характеристикой оператора интегрирования дробного порядка М > 0 по координате состояния, т.е. бЕ^ ] = Р~м(\,1).
2. Операторы дифференцирования. Будем рассматривать оператор дифференцирования по времени с учетом значения функции в начальный момент [6-10]:
Б0^к(в, х) = ^в1 + ¿(вЖв, х), (15)
а также операторы дифференцирования первого и второго порядков по координате состояния:
ВАКв, х) = , Е2, И(в, х) =д-Щхх1 (16)
а+|х ( , ) дх а+|х ( , ) дх2 v 1
где И(в,х) е Ь2(О), Q = Тхй, Т = [0,Г], й = [а,Ь] с Я .
Спектральные характеристики операторов дифференцирования и Еа+^х,
определенные относительно базисной системы {е(Ч0,ч,в,х)}°^=0 пространства Х2(О), вычисляются по определению спектральных характеристик линейных операторов и обозначаются Р(2,2) =(Р1Мк ), 3!(2,2) = (%11М1К2 ) соответственно:
Ры К К =( е&, Ь,в, х), де°0±в, х) 1 +{е(ч, Чв х), еу (л, Кв х) \ (в),
V дв )ь1(О)
3,0,= [еЦ^Ч,в,х) Ш-/0,Квх) 1 , ^ч^Л = 0,и....
\ дх )ьг(О)
Спектральная характеристика оператора дифференцирования Е2+|х обозначается
через Р (2,2) и определяется по свойству спектрального преобразования композиции линейных операторов:
3 1(2, 2) = б[Е(2+| х ]= Б[е«+| х ° Еа+\ х ]= Б[е«+|х ] х ]=31 (2,2) ' 31 ^2) .
Если функции базисной системы пространства (О) представляются в виде {е(Ч0,ч,в,х) = <(Ч0,в)• р(\,х)}°°г. =0, т.е. в виде произведения функций базисных систем пространств ¿2 (Т) и ¿2 (й), то
Р(2,2) = Р(1,1) ® £(1,1), ^(2,2) = Е (1,1) ®ЗД,1), Зп(2,2) = £(1,1) ®3П(1,1). (17)
В приведенных выражениях Р(1,1) - спектральная характеристика оператора дифференцирования по времени с учетом значения функции в начальный момент, определенная относительно базисной системы {з(',в)}^0 пространства ¿2 (Т), а ^ (2,2) и Зп (1,1) - спектральные характеристики операторов дифференцирования первого и второго порядка по координате состояния, определенные относительно базиса {р(\,х))}^=0 пространства ¿2(О) [5, 9, 10].
Композиция операторов дифференцирования (15)-(16) и дробного интегрирования (10)-(11) порождает операторы дробного дифференцирования [1, 5].
Приведем здесь только два свойства спектрального преобразования операторов дробного дифференцирования.
а) Свойство спектрального преобразования операторов дробного дифференцирования (модифицированного) у -го порядка (0 <у< 1) по времени при заданных начальных условиях:
£
Б'вМв, х)
_л\в в х)=к_л\в0( х)
= Ру (2,2) • Н (2,0) _ ^(1,0; 00) ® Н^, (1,0)
(18)
' Квкв, х)
к(в0,х )=к0( х )
= Р\ (2,2) • Н (2,0) _ (р-(1-Л) (1,1) • р(1,0; в,))® Н, (1,0)
где Ру (2,2) = Ру (1,1) ® £(1,1) ( Р\ (2,2) = Р\ (1,1) ® £(1,1) ) - спектральная характеристика оператора (модифицированного оператора) дробного дифференцирования первого порядка по времени, определенная относительно базисной системы разложения {з('0 ,в)}°
пространства ¿2 (Т), а Н ^0(1,0) (Н0 (1,0)) - спектральная характеристика начального
условия (х) (к0 (х)), определенная относительно базисной системы разложения
{р(\, х)}^0 пространства ¿2 (О) .
Ь) Свойство спектрального преобразования операторов дробного дифференцирования (модифицированного) а -го порядка (1 <а< 2) по координате состояния при заданных краевых условиях первого рода:
8
5
Ба И(9, х)
К-а\х (в,а)=кг_а,а (в), й2_,
~КЪ)=кг_а> „ (в)
= 3аи(2,2) • Н (2,0)-3,(2, 2) ■
2-а Ъ
(1,0) ® р(10; Ъ) - 1^2-а а (1,0) ® р(1,0; а)] (19)
(
5
* Б" И(в, х)
Л
,а)=Иа (в), Ь(в,Ъ)=\(в)
= 3,1* (2,2) • Н (20) - 3" (2,2) • [щ (1,0) ® р(1,0; Ъ) - [ (1,0) ® р(1,0; а)] где 3а1(2,2) = Е ®(- 3г(1,1)3[ (1,1)РГ(2-а)(1,1)), 3^. (2,2) = Е ®(- Р1-(2-а)(1,1)31(1,1)3Г (1,1)), 3 (2,2) = Е ®3 (1,1), 3 (2,2) = Е ® р-(2-а) (1,1)3 (1,1) - спектральные характеристики операторов (модифицированных операторов) дробного дифференцирования второго и первого порядков по координате состояния, определенные относительно базисной
системы разложения {р(ч,х)}"=0 пространства Х2(О), а Щ2-а\ъ(1,0) и щ2_а|а(1,0) (щ(1,0) и [(1,0)) - спектральные характеристики краевых условий ь(в) и Л2_а|а(в)(къ(в) и ка (в)), определенные относительно базисной системы разложения {з(/0 ,в)}^ пространства (Г).
3. Гиперматричные передаточные функции ОРП и СРП и переходные блоки
В дальнейшем, по мере необходимости, при описании систем управления в спектральной области, спектральные характеристики любых линейных операторов (5) будем называть гиперматричными передаточными функциями (ГМПФ).
Как известно [22,23] одномерная ИПФ О(в,т0, х,£0) ОРП является решением краевой задачи
а (в, х)
дги(в,х)
дв2
+а (в, х)
ди(в,х) дв
= Ъ (в, х)
д 2и(в, х)
дх2
+ Ъ2(в, х)
ди(в, х) дх
+
+ Ъ (в, х)и(в, х) + ё(в-т0)5(х -£0); а < х < Ъ, в> 0;
и(0, х) =
ди(0, х) дв
= 0, х е [а, Ъ];
(20)
(21)
а(в, а)и(в, а) + (5(0, а) ди(в, ^ = 0, в> 0;
дх
а(в, Ъ)и(в, Ъ) + ((в, Ъ) дивЪ = 0, в> 0
дх
(22)
(23)
и существенно зависит от краевых условий (22)-(23). Преобразование задачи (20)-(23) в спектральную область позволяет найти ГМПФ ОРП. Кроме того, формулы связи для ДНПФ последовательного, параллельного соединения и соединения с обратной связью в спектральной области для систем с сосредоточенными параметрами [6,7] остаются
справедливыми и для ГМПФ последовательного, параллельного соединения и соединения с обратной связью систем с распределенными параметрами [11].
Пусть регулирование распределенной системы [22, 23] осуществляется с помощью сосредоточенного регулятора, который измеряет состояние распределенной системы в одной точке одномерной области определения выходного сигнала и производит регулирующее воздействие на распределенный блок также в одной точке области определения его входного сигнала.
Формирование математической модели такой системы управления в спектральной области связано с введением двух переходных блоков, у которых пространственная размерность входного сигнала не совпадает с пространственной размерностью выходного сигнала. Для одномерного случая эти переходные блоки описываются ИПФ 8(х — а) и 8(а—^), а их ГМПФ имеют вид
Жа (2,2) = Е ® Б а (24)
и
где
(2,2) = Е ® Б, (25)
ВД, г) = \\8(а—^)р(], х)р(/,£У£аХ = { р(], х)Л • р(г,а), (26)
рр □□ □
г) = \\8(х — Р)р(],х)р(г,^¿х = р(], Д)\р(г,. (27)
рр * * *
У
рр
□ □ □
Заметим, что более общие модели переходных блоков рассмотрены в работах [22, 23]. Для них также легко находятся их ГМПФ.
3. Моделирование диффузионных процессов с применением пакета расширения МЬ8У_8М+Ма1Иеаа
Пример 1. Моделирование процесса аномальной диффузии Будем решать спектральным методом краевую задачу для уравнения аномальной диффузии [21]
Бо1вк(в, х) = а( х)Баа+1 к(в, х) + Е (в, х), (28)
которое получается из классического дифференциального уравнения диффузии
^ = а( х) ^ + 8 (в, х)
дв дх
заменой производной второго порядка по пространственной координате на дробную производную Римана-Лиувилля порядка 1 <а< 2 и производной по времени первого порядка на дробную производную Римана-Лиувилля порядка 0 < у < 1.
Здесь, как и в стандартном уравнении диффузии к(в, х) - функция концентрации вещества на отрезке а < х < Ъ; % (в, х) - функция источника/стока на данном отрезке; а(х) - коэффициент диффузии. Из физических соображений следуют условия 1 <а< 2, 0 < у < 1 и а(х) > 0. Также имеем \ (в,х) = (0,х) = /(х) для а < х < Ъ и граничные
у I в=0 у
условия первого рода К"0, х)| х=а = К"0, а) = /а (в) и К2-а(0, х)| х=Ъ = К2-а(0, Ъ) = Л (в) для всех в> 0, где ЬЩв(в, х) = ¡^к(в, х) , Ь^ (в0, х) = ¡^Кв, х) [5].
В качестве примера рассмотрим начально-краевую задачу, которая заданна в виде:
Б1;и(в, х) = а( х) Б0+7 И(в, х) + % (в, х),
К-У (0, х) = 0, к2-а(в, а) = 0,
0 < х < Ъ, 0 <в< 1,
(29)
где % (в, х) =
И2_а (в, Ъ) = 0745 •в2, 0 <в< 1,
Г(3) 2 „2, .2 „^ _П13),7
-в13х2 -в2х2 и а(х) = ■
-х
Г(2.3) Г(3)
Точным решением начально-краевой задачи (29) является функция к(в, х) =в2 х2, что может быть проверено ее прямой подстановкой в (29) и использованием формулы дробного дифференцирования [1, 2, 5].
Будем решать эту задачу спектральным методом. Для этого преобразуем заданную краевую задачу в спектральную область. Для спектрального преобразования функций времени и состояния в качестве базиса пространства , где Q = Г х О, Г = [0, ^],
Q = [0, Ъ], выберем такую ортонормированную систему
{е(}0, гг ,0, х) = ф(}0, в) • р(К, , (30)
что системы функций {з(/0 ,в)}}, {р(^, х)}} являются базисными системами пространств
Х2 (Г), (О) соответственно. Применим спектральное преобразование к начально-краевой задачи (29). Тогда
5
Б0°вК(в, х)
= 5
V 30 (0,х)=0
а( х) • В 7 К(в, х)
К3\х(в,Ъ)=0.7450 ,Ь033\х(в,0)=0
+ 5 [# (в, х)]
(31)
Рассмотрим левую часть равенства (31). Используем свойство спектрального преобразования операторов дифференцирования по времени при заданных начальных условиях (18), левую часть (31) представим в виде:
5
Бвчв, х)
= Р07 (2,2) Н (2,0), (32)
^ Ао.з| в (0,х)=0
где Р0.7(2,2) спектральная характеристика оператора дробного дифференцирования по времени Римана-Лиувилля с учетом значения функции в начальный момент, определенная относительно базисной системы (30), для которой, согласно (17), справедливо представление
Р0 . 7 (2,2) = Р0 7 (1,1) ® Е(1,1), где Р07(1,1) = Р(1,1)Р 03(1,1) - спектральные характеристики оператора дробного дифференцирования по времени Римана-Лиувилля с учетом значения функции в начальный момент времени, определенная относительно базисной системы {з(г0 ,в)}}3;
Е(1,1) - двумерная единичная матрица. Через Н(2,0) обозначена спектральная характеристика функции к(в, х), определенная относительно той же базисной системы
(30).
Рассмотрим теперь правую часть равенства (31). Используем свойство спектрального преобразования операторов дробного дифференцирования по состоянию при заданных краевых условиях (19) и свойство композиции операторов (7), правую часть
(31) представим в виде:
5
'и в, х
к03, х (в,Ь)=0.745в2,й03, х (в,0)=0
а( х) • Б'7 И(в, х) + С (2,2)^ (1,0) ® р(1,0; Ь)) + 0(2,2)
+ 5 (в, х)] = 5(2,2) Н (2,0) +
(33)
Ьу
В (33) 5(2,2) = Е(1,1) ® (— А(1,1) • 3! (1,1) • 3[ (1,1) • Р—°-3 (1,1)) и
С(2,2) = Е(1,1) ®(Л(1,1) •31(1,1)) - четырехмерные гиперквадратные матрицы, где А(1,1) -спектральная характеристика оператора умножения, определенная относительно базисной системы {р(^, х)}^; 3 (1,1) - спектральная характеристика оператора дифференцирования
по координате состояния, определенная относительно базисной системы {р(\, х)}°; Е(1,1) - двумерная единичная матрица. Н (2,0) - спектральная характеристика функции к(в, х), определенная относительно базисной системы {з(/0 ,в)- р(г, х)}^ . ^ (1,0) - спектральная
2
характеристика краевого условия (Ь, х) = УЬ (х) = 0 745в , определенная
относительно базисной системы {з(/0,в)}°. Через р(1,0;Ь) обозначена матрица-столбец функций базисной системы {р^, х)}° в точке х = Ь.
С учетом введенных обозначений уравнение (31) можно переписать следующим образом
[р07 (2,2) - Б(2;2)Н (20) = С(2,2) + С(2,2)(% (1,0) 0 р(1,0; Ъ)). (34)
Выражая спектральную характеристику Н(2,0) из уравнения (34), получаем
Н (20) = [Р01(2;2) - Б(2,2)]"1{а(2,2) + С(2,2)(Х (1,0) 0 р(1,0\ Ъ))}. (35)
Для получения решения задачи (29) в пространстве функции времени и состояния требуется применить формулу обращения:
Ь1-\ Ь2 -1
к(в, х) = £ £ Н„-№,0) • р(/, х), (36)
]=0 ¿=0
где - координаты спектральной характеристики Н(2,0) , 0 е [0, ?], х е [0, Ъ].
Для решения этой задачи используем пакет расширения 8рек1х_8М СКМ МаШсаё [5, 12]. В качестве базисных систем {р(/0,в)}}, {р(\,х)}} будем использовать полиномы
Лежандра.
Решение задачи:
Листинг 1.
N1 := 4 <-порядок усечения матрицы ДНПФ по времени
N2 := 4 <-порядок усечения матрицы ДНПФ по состоянию
1 := 1 Ь := 1 Ы := 20 Ь2 := 10
1^(4) := 0.745 ■ Ь <- краевое условие
а(х) := —:—- ■ х Г(3)
е1(е) :=
(Ш-.е13
\ Г(2.3)
- в2 §2(х) := х2
02 := ЭКХРР1(^,Ю,Ь) 01 := ЗКХРР1(^1 ,N1 Д) т
РМ := 01 ■ 02 <-ДНСХ входного распределенного воздействия. 0:=РС2(та) А := ЭУ2РР1(а,К1 Д) А := ЭУ2РР1(а,N2,Ь)
:= Э1рРР1(М1 ДДЗ) 1р2 := 31рРР1(Ю,Ь ДЗГ) Ч^ := ЗИХРР1|Т4/1 ,N1 Д| <-НСХкраевого условия.
N51 := ЭКВРР1(]_1,М1 Д)
ЫВх := ЗМВРР1('Ь2,К2,Ь') Е1 := ¿с1епШу(М1) Е2 := ¿с1епШу(Н2Г)
<- БС полиномов Лежандра.
Р1 := 5Р1РРЗ(Ь,М2)
II := 5ПРР1(^т) ТЛТ := 5Р1РР1(1.Х1} ■ 1р1 ШХ
< ДНПФ эпеметарньп -ЗЕеньеЕ.
:= А ■ I
Р1 РГ
чн)
1 := О- XI - 1 := ХВ^1 ] := О.. N2 - 1 А^ 0 := КВ*^
:= (ТР(^Т:Е2) + ГР(Е1Л\1Х))-1
<-ГМПФДУ. Здесь 1Р(А_В)- программа вычисления «якорного произведения матриц А и В. т
К¥ := - АЬ <-матрицы ДНСХ краетаи значений.
ОК := [ТР[Е1ГА ■ (Р1 ■ 1}] ■ ГС2(КТ)]
<- ГСМ ДНСХ краевых значений. Здесь ГС2(А) - программа вычисления птерстолбцовой матрицы по матрице А.
Н-5^ := № ■ (ОК + О) <- ГСМ ДНСХ выходного распределенного сигнала.
Н := <- матрица ДНСХ выходного распределенного сигнала.
Коней листинга 1.
Заметим, что если в начально-краевой задачи (29) вместо дробных производных Римана-Лиувилля использовать модифицированные производные Римана-Лиувилля [1, 2, 5], то она примет вид:
*О(0, X) - а( х)-Ж хОД х) = % (0, х),
(37)
(38)
й(0, х) = 0, 0 < х < Ъ, И(в,а) = 0, 0 <0< г, И(в,Ъ) = в\ 0<0<г,
где % (0, х) = -Г3)- 2-^х 2 и а(х) = Г13) х17. Г(2.3) Г(3)
Применяя спектральное преобразование к начально-краевой задачи (37), получим
[е 7(2,2) - А(2,2р;; (2,2)]н (2,0) = С(2,0) +
+ ¿(2,2)317 (2,2)(^ъ (1,0) ® р(1,0, Ъ)). .
Выражая спектральную характеристику Н (2,0) из уравнения (38), найдем
Н (2,0) = [Р,07 (2,2) - ¿(2,2)31; (2,2)]-1 • {а(2,0) + + ¿(2,2)31.7 (2,2)(^ъ (1,0) ® р(1,0, Ъ))}. В (39) 3Ц (2,2) = £(1,1) ®(-3.7(1,1) -3[ (1,1))= Е (1,1) ®(- Р1-0-3(1,1)31(1,1)-3 (1,1)) и 3 (2,2) = £(1,1) ®3[ (1,1) - четырехмерные гиперквадратные матрицы, где 3 (1,1) -спектральная характеристика оператора дифференцирования по координате состояния, определенная относительно базисной системы {р(\, х)}°, а £(1,1) - двумерная единичная
матрица. Н(2,0) - спектральная характеристика функции к(в, х), определенная относительно базисной системы {(з(/0,0)- р(\, х)}°. . ^ (1,0) - спектральная
(39)
2
характеристика краевого условия первого рода ^(0,Ъ) = УЪ (0) = 0 , определенная
относительно базисной системы {p(i0,в)}^. Через p(1,0; b) обозначена матрица-столбц функций базисной системы {p(i, x)}°° в точке x = b .
Точным решением этой задачи так же является функция И(в, x) = в2x2, которая является решением задачи (29). Результаты решения задачи (37) спектральным методом совпадают с результатами, приведенными в листинге 1. Это можно показать, используя формулу связи [5] между операторами дробного дифференцирования Римана-Лиувилля и модифицированными операторами дробного дифференцирования Римана-Лиувилля. Для
нашей задачи имеем: Dlgh(e,x) = —h(°'x)—+"Dlgh(e,x), 0 < у < 1;
1 Г (1-у) •в 1
Da. h(0, x) =-+-x 4 w , +*D", h(d,x), 1 < a < 2.
°+lx v ' 7 гп.„л. va va-1 0+lx v ' ''
Н(в, х) к(1\в,х) ПуОух* Г (2-а) • ха
Пример 2. Моделирование нейтронной кинетики ядерного реактора в одногрупповом диффузионном приближении
Рассмотрим описание переноса нейтронов в одногрупповом диффузионном приближении, относящегося к случаю, когда коэффициент диффузии можно принять не зависящим от пространственных переменных [19, 21]:
1 дФ^т) = т 2ф(0, г) _ £ а тще, г)+ К де (40)
M
+(1 - ^эффК (в, - (в, те, Г)+Z ЪС (в, П
i=1
В уравнении (40): Ф - плотность потока нейтронов; С1 - концентрация ядер-предшественников запаздывающих нейтронов I -той группы; V - средняя скорость нейтронов; Б - коэффициент диффузии; Еа - сечение поглощения; Д,фф - эффективная
доля запаздывающих нейтронов; - постоянная распада ядер-предшественников / -ой группы; = vZ/ /Еа - коэффициент размножения нейтронов для бесконечной среды, где Е/ - сечение деления, а V - среднее число нейтронов.
Перепишем это уравнение для плотности нейтронов N (с учетом того, что
Ф = N •V, V - скорость нейтрона) совместно с уравнением для концентрации предшественников запаздывающих нейтронов. Получим:
SN LL Гр-Дфф_+_ LB2
дв l-кэфф • (1 + L2B2)
-V2 N +
V l l• кэфф • (1 + L2B2) j
N + , (41)
i=1
SC вэфф
дС =-Л,С N. (42)
де I
В уравнениях (41) и (42): N - плотность тепловых нейтронов (мощность); } = Б / Е - квадрат длины диффузии тепловых нейтронов; В2 - материальный
лапласиан; Д
эфф
эффективная доля запаздывающих нейтронов i -ой группы;
кэфф = кш /(1 + Ь2Б2) - эффективный коэффициент размножения нейтронов; l = 10 /k,
О ' "'эфф
среднее время жизни нейтронов, где 10 = 1Ш /(1 + Ь2Б2) - время жизни нейтронов в среде конечных размеров, а 4, = 1/У - среднее время жизни нейтрона в бесконечной среде;
Р = (kэфф "1)/ kэфф - реактивность.
Линеаризуем диффузионные уравнения (41), (42). Для этого представим значения
каждой из функций N(в, г), С(в, г) в отклонениях от соответствующих величин в стационарном состоянии, т.е.
N (в, г) = N (О, г) + (в, г) = N + ¿V; С, (в, г) = Сг (О, г) + (в, г) = Сг0 +
и учтем, что в стационарном состоянии V2N(r) + B2N(r) = О. Тогда для плоского реактора (одномерный случай) в отклонениях от стационарного состояния их можно записать, с учетом начальных и краевых условий, в виде
д~(в, х)
дв
I - (1 + Ь В О
д 2~(в, х)
дх2
+ В2 ~(в, х)
+ -
Д
эфф
I
~(в, х) -
-11 Дэфф~ (в, х) р(в, х);
I
г=1
I - N°
(43)
д~ (в х) ДэффУтах дС (в х) + (в, х) = Д-^ ~ (в, х);
дв г I - Стах
О
г = О,1,... 6; О < х < Ъ;в> О;
~ (0, х) = ~о (х) = 0, ~ (0, х) = ~О (х) =О;
~(в,О) = ~(в, Ъ) = О,
где ~ (в, х) = ¿У(в, х) / У0т1х - безразмерная нормированная плотность потока нейтронов, ~(в, х) = ЙС(в, х)/ Стах - нормированное отклонение концентрации ядер-
предшественников запаздывающих нейтронов г - ой группы.
Решение начально-краевой задачи (43) выполним с использованием спектральной формы математического описания [9-12]. Для спектрального преобразования функций времени и состояния в качестве базиса пространства Ь2 (<2), где Q = Т х О, Т = [0, ^ ],
0 = [О,Ъ], выберем такую ортонормированную систему {е(/0,^,в,х) = ^(/0,в)-р^,х)}^ 0, что системы функций {з(/0 ,в)}°=0 и {р(\, х)}°=0 являются базисными системами
пространств Ь2 (Т) и Ь2 (О) соответственно.
Применим спектральное преобразование к системе уравнений (43). Так как ГСМ НСХ частных производных, входящих в систему уравнений (43), с учетом начальных и
д~(в, х)"
краевых условий,
вычисляются по правилам:
£
дв
= Р(2,2) N (2,0),
£
£
дсг(в, х)'
дв _
" М0(х)
= Р(2,2)С, (2,0), £
д2 ~(в, х)
дх2
= 3П(2,2) N(2,0),
N
Ф(в,х)
= А(2,2)0(2,0),
где
а ГСМ НСХ
Р(2,2) = Р(1,1) ® Б(\,\),
3(2,2) = Е(1,1) ® (-3(1,1) 3 (1,1)), А(2,2) = Б(1,1) ®А(1,1), то краевая задача (43) в спектральной области записывается в виде
Р(2,2) N(2,0) -
Ь2
[3! (2,2) N(2,0) + В2 (2,0)]+ ^^ #(2,0) -
I • (1 + Ь2 В2)
1 6 ~ 1
- - ХДЭфФССг (2,0) = - А(2,2)0(2,0);
1 1=1 1
(44)
тах
Р(2,2)Сг (2,0) + ЛС (2,0) = Р\ • N(2,0).
I • С
тах 0
Решая систему матричных уравнений (44), находим ГСМ НСХ безразмерной нормированной мощность реактора
N(2,0) = (2,2)0(2,0), (45)
где
(2,2) =1 | Р(2,2) -
Ь2
+
I | I • (1 + Ь2В2)
^Л^^ФФ[Р(2,2) + АгЕ(2,2)]-1 I А(2,2)
[з ,(2,2) + В2 Е ]+^эфФ Е +
Жтах 6
0 V1 1 ДЭ
(46)
I • С" г=1
0 г 1
- ГМПФ кинетики нейтронов ядерного реактора, найденная для модели с шестью группами запаздывающих нейтронов, а
0(2,0) = £[ф(в, х)] (47)
- ГСМ НСХ входного воздействия.
Для получения решения задачи (43) в пространстве функций времени и состояния применим формулу обращения
Ь -1 ь2 -1
u(в, х) = £ £ Ы],г ' , в) ' Р(1, х) ,
]=0 г=0
где и - координаты ГСМ НСХ и(2,0), в е [0, г], х е[0, Ь].
(48)
Для решения этой задачи используем пакет расширения 8рек1х_8М СКМ МаШеаё [5, 12]. В качестве базисных систем {р(/0,в)}^, {р(г,х)}^ будем использовать полиномы Лежандра. Листинг программы расчета здесь не приводится.
<
Результаты моделирования по рассмотренным алгоритмам, заданной реактивности pß) = 0.1 -ß -1ß- 40), а также параметрам, которые определяют кинетику ядерного реактора [19, 21]:
ж
b = 850, L = 2.23, B =-, l = 0.001, ß = 0.0065,1 = 0.0124, 1 = 0.0305,
2-b 4 2
1 = 0.111,1= 0.301,1= 1.14,1 = 3.01,ß = 0.0002145, ß = 0.0014235, ß = 0.0012740, ß = 0.0025675, ß = 0.0007475, ß = 0.0002730,
приведены на рис.1. На рис.1. а) показаны сечения поверхности безразмерной нормированной плотности потока нейтронов в плоском реакторе во времени, а на рис. 1. б) показана сама поверхность относительной плотности распределения потока нейтронов во времени и пространстве.
2.3 -
1.3
hi.3 ^ 1
4+ 0.3
Ks
0'5 0 60 120 180 240 300
Рис. 1
4. Моделирование системы управления ядерной энергетической установкой с распределенными параметрами с применением пакета расширения МЬ8У_8М+МаШсаа
Основной задачей автоматического управления ядерной энергетической установкой является управление нейтронной мощностью. Применение спектрального метода для математического моделирования САУ мощностью ядерного реактора, уравнения кинетики которого описываются точечной моделью, рассмотрено в работе [19]. Математическую модель САУ ядерной энергетической установкой, рассмотренной в работе [19], можно модифицировать под реактор, который описывается уравнениями кинетики с распределенными параметрами. Спектральная схема такой модифицированной САУ ядерного реактора приведена на рис 2.
Рис. 2
В этой модели кинетика нейтронов ядерного реактора описывается ГМПФ (46), которая найдена для распределенной модели с шестью группами запаздывающих нейтронов.
Тепловые процессы, управление мощностью реактора, а также дробный ПИД регулятор и эффекты запаздывания учтены ГМПФ: Щ (2,2) = Щ (1,1) 0 Е, Жрс(2,2) = Жрс(1,1) 0Е, Жпид(2,2) = Жпид(1,1) 0Е, Жсо(2,2) = Жсо(1,1) 0 Е соответственно.
Эти ГМПФ сами выражаются через ДНПФ Щ (1,1) = Щ (г, г) = -а- Т0Л -[тосР(г, г) + Е]-1,
Щрс(1,1) = Щрс(г, г) = кпр[тпрР(г, г) + е] - Р~\г, г),
Щщщ (1,1) = Щщд (г, г) = кхЕ+£2За(г, г) + к}Р~м(г, г), (1,1) = т~00(г, г), рассмотренные в работе [19].
Кроме того спектральная схема САУ ядерного реактора содержит два переходных
блока Ща(2,2) и Щр(2,2) (см. формулы (24)-(27)), в которых параметр а = ( = Ь/2.
По спектральной схеме САУ ядерного реактора, ГМПФ звеньев, входящих в ее состав, находим четыре ГМПФ:
(2,2) = [Е 0 Е + Щ (2,2) - Щрс(2,2) - Щшд (2,2) - т"^ (2,2)Щ (2,2)}'
хЩ(2,2)-Щс(2,2)-Щпид (2,2),
(2,2) = [Е 0 Е + Жкин (2,2)Щр (2,2)Щ (2,2) + + Щро (2,2) - Жпид (2,2) - т"00 (2,2)Ща (2,2]* (2,2Щ (2,2), (2,2) = [Е 0 Е Щ (2,2) - Щрс (2,2) - ЖпЩ (2,2) - т"00 (2,2): Щ(2,2) - Щин(2,2)У - Щ(2,2) - Щрс(2,2) - Щид (2,2),
х
:
:
(49)
(50)
(51)
(2,2) = [Е 0 Е Щ (2,2) - Щрс (2,2) - Щид (2,2) - т"00 (2,2) хЩ((2,2)-Щкин(2,2)У Щ(2,2),
:
где
^(2,2) = [к ® Е + Жшн (2,2) Жд (2,2) Ж (2,2)^(2,2)]^ (2,2)^(2,2),
Ж (2,2) = [Е ® Е + Жд (2,2)Жт (2,2) • Ж (2,2) • Жшн (2,2)]^ • ЖД (2,2) .
По ГМПФ (49) и (50), ГСМ НСХ приращения мощности потока нейтронов ~ и
реактивности ртзм, вносимой внешним возмущающим воздействием, находим ГСМ НСХ для относительного приращения мощности потока нейтронов п
N(2,0) = (2,2) • ¿/(2,0) + (2,2) • £(2,0) . (53)
По ГМПФ (51) и (52), ГСМ НСХ приращения мощности потока нейтронов ~ и реактивности р , вносимой внешним возмущающим воздействием, находим ГСМ НСХ для относительного приращения реактивности р
(7(2,0) = ЖР (2,2) • ¿7(2,0) + ЖрР_ (2,2) • £(2,0) . (54)
Для решения этой задачи используем пакет расширения Spektr_SM СКМ Mathcad [5, 12]. В качестве базисных систем \р(}0,0)}}, \р(Ч,Х)}Г будем использовать полиномы
Лежандра. Листинг программы расчета здесь не приводится.
Некоторые результаты моделирования САУ ядерной энергетической установки с распределенными параметрами спектральным методом (см. раздел 3, пример 2) по рассмотренным алгоритмам и заданных управлении ~(0) = 0.05 •(1(0-20) - 0.9 • 1(0-75)) и возмущении рвозм (0) = 0.1 • Д • (1(0 - 40)), а также параметрах САУ ядерной энергетической
установки, заданных в [19], приведены на рис.3. На рис. 3. а) показаны сечения поверхности безразмерной нормированной плотности потока нейтронов в плоском реакторе во времени, а на рис. 3. б) показана сама поверхность относительной плотности распределения потока нейтронов во времени и пространстве.
0.07
0 055 0 60 120 180 240 300
Рис. 3
В заключение заметим, что в данной статье спектральный метод применяется для моделирования распределенных и дробно-распределенных процессов и систем управления. Для моделирования таких процессов и систем модифицирован и применен пакет расширения MLSY_SM СКМ Mathcad. Этот пакет, а также пакеты расширения MLSY_SM, Spektr_SM+Simulink+Matlab, Spektr_SM+VisSim+Mathcad СКМ [12-19] позволяют проводить не только детерминированный, но и стохастический анализ распределенных САУ различных классов, а также проводить параметрическую оптимизацию параметров распределенных САУ. В частности можно исследовать влияние степени дробности операторов дробного интегрирования и дробного дифференцирования на переходные процессы, протекающие в САУ [5] и по критериям оптимальности решать задачу оптимизации параметров регулятора.
Библиографический список
1. Самко С.Г., Килбас А.А., Маричев О.И. Интегралы и производные дробного порядка и некоторые их приложения. Минск: Наука и техника. - 1987. - 688 с.
2. Васильев В.В., Симак Л.А. Дробное исчисление и аппраксимационные методы в моделировании динамических систем. - К.: НАН Украины, 2008. - 255с.
3. Потапов А.А., Гильмутдинов А.Х., Ушаков П.А. Фрактальные элементы и радиосистемы: Физические аспекты / Под ред. А.А. Потапова (Библиотека журнала «Нелинейный мир»: Научная серия «Фракталы. Хаос. Вероятность»).- М.: Радиотехника, 2009. - 200 с.
4. Бекмачев Д.А., Потапов А.А., Ушаков П.А. Проектирование фрактальных пропорционально-интегрально-дифференциальных регуляторов дробного порядка // Успехи современной радиоэлектроники. 2011. №5. С. 13 - 20.
5. Рыбин В.В. Моделирование дробных нестационарных систем управления в СКМ спектральным методом // Вестник Московского авиационного института. - 2011. Т. 18. № 6. - С.102 - 118.
6. Солодовников В.В. и др. Расчет систем управления на ЦВМ: спектральный и интерполяционный методы. - М.: Машиностроение, 1979.- 664 с.
7. Солодовников В. В., Семенов В. В. Спектральная теория нестационарных систем управления. М., Наука, 1974.- 336 с.
8. Нестационарные системы автоматического управления: анализ, синтез и оптимизация / Под ред. К.А. Пупкова и Н.Д. Егупова. - М.: Издательство МГТУ им. Н.Э. Баумана, 2007. - 632с.
9. Пантелеев А.В., Рыбаков К.А., Сотскова И.Л. Спектральный метод анализа нелинейных стохастических систем управления - М.: Вузовская книга, 2006. - 392 с.
10. Пантелеев А.В., Рыбаков К.А. Прикладной вероятностный анализ нелинейных систем управления спектральным методом. - М.: Изд-во МАИ-Принт, 2010. - 160 с.
11. Клевцов Ю. А. Спектральный метод моделирования и идентификации объектов с распределенными параметрами. Диссертация на соискание ученой степени кандидата технических наук. - Киевский политехнический институт. Киев,1984. - 97 с.
12. Рыбин В.В. Моделирование нестационарных непрерывно-дискретных систем управления спектральным методом в системах компьютерной математики.. - М.: Изд-во МАИ, 2011. - 220 с.
13. Рыбин В.В. Разработка пакета расширения MLSY_SM СКМ Mathcad в биортогональных вейвлет-базисах // Электронный журнал "Труды МАИ" - 2009, № 33. - http://www.mai.ru
14. Рыбин В.В. Описание сигналов и линейных нестационарных непрерывных систем управления в базисах вейвлетов. // Электронный журнал "Труды МАИ" - 2003, № 10. -http://www.mai.ru
15. Рыбин В.В. Разработка и применение пакетов расширения MLSY_SM СКМ Mathcad, Maple, Mathematica, Matlab.// Электронный журнал "Труды МАИ"- 2003, № 13. -http://www.mai.ru
16. Рыбин В.В. Разработка и применение пакета расширения Spektr_SM СКМ Matlab.// Электронный журнал "Труды МАИ" - 2003, № 13. - http://www.mai.ru
17. Рыбин В.В. Разработка пакета расширения МЬ8У_БМ СКМ МаШсаё в базисах Добеши М-го порядка // Электронный журнал "Труды МАИ"- 2009, № 33. - http://www.mai.ru
18. Рыбин В.В. Разработка пакета расширения МЬ8У_БМ СКМ МаШсаё в проекционно-сеточных финитных базисах // Электронный журнал "Труды МАИ" - 2010, № 41. -http://www.mai.ru
19. Рыбин В.В. Моделирование САУ ядерной энергетической установкой в СКМ спектральным методом // Электронный журнал "Труды МАИ" - 2012, № 50. -http://www.mai.ru
20. Петухов А.А., Ревизников Д.Л. Алгоритмы численного решения дробно-дифференциальных уравнений // Вестник Московского авиационного института. -2009. Т. 16. №6. С.228 - 234.
21. Халимончук В.А. Динамика реактора с распределенными параметрами в исследованиях переходных режимов эксплуатации ВВЭР и РБМК. - К.: Основа, 2008. - 228 с.
22. Бутковский А.Г. Характеристики систем с распределенными параметрами (справочное пособие). - М.:Наука,1979. - 224 с.
23. Рапопорт Э.Я. Структурное моделирование объектов и систем управления с распределенными параметрами. - М.: Высшая школа, 2003 г. - 299 с.
Сведения об авторах
Рыбин Владимир Васильевич, доцент Московского авиационного института
(национального исследовательского университета), тел.: +7 499 158-48-11, е-таП: уу-
ribin@mail.ru
