Моделирование теплотехнологических объектов в классе дробно-иррациональных передаточных функций Текст научной статьи по специальности «Математика»

УДК 62-523.8

МОДЕЛИРОВАНИЕ ТЕПЛОТЕХНОЛОГИЧЕСКИХ ОБЪЕКТОВ В КЛАССЕ ДРОБНО-ИРРАЦИОНАЛЬНЫХ ПЕРЕДАТОЧНЫХ ФУНКЦИЙ

Е.Б. КАРИКОВ1 В.В. МИШУНИН2 В.Г. РУБАНОВ1 Ю.А. ГОЛЬЦОВ1

В статье описаны методы идентификации систем с распределенными параметрами, использующие частотные и временные характеристики. Для идентификации по временным характеристикам был модифицирован метод идентификации Симою. Приведены результаты идентификации различными методами на примере решения задачи нагрева полубесконечного тела.

Ключевые слова: идентификация, дробно-иррациональная передаточная функция, метод наименьших квадратов, Симою.

e-mail: rubanov@intbel.ru

Качество процесса управления напрямую зависит от адекватности математической модели, используемой при синтезе алгоритмов управления. Многие из теплотехнологических объектов управления относятся к классу объектов с распределенными параметрами и описываются дифференциальными уравнениями в частных производных. Для использования классических методов синтеза закона управления приходится прибегать к аппроксимации данной модели стационарными дифференциальными уравнениями с постоянными коэффициентами, что не позволяет учесть всех особенностей динамики процессов, протекающих в системах с распределенными параметрами. Повысить качество конечного продукта и снизить затраты энергии позволяет использование системы управления, построенной на основе математического описания, отражающего физическую сущность тепловых процессов, описываемых дифференциальными уравнениями в частных производных. Аппроксимация таких уравнений в области комплексного переменного передато ч-ными функциями целого порядка приводит к модели с невысоким уровнем адекватности. В то время как аппроксимация объекта управления дробно-иррациональными передаточными функциями наиболее полно отражает распределенный характер процессов, протекающих в объекте управления, что позволяет достичь более высокой степени адекватности, а, следовательно, и то ч-ности управления процессом.

Получение математической модели объекта автоматизации с дробным порядком может осуществляться различными методами идентификации, например, по частотным характеристикам или по кривым разгона, снятым экспериментально, при этом точность описания будет зависеть от выбора метода идентификации.

Для использования метода идентификации по частотным характеристикам [1] передаточную функцию объекта можно записать в виде отношения полиномов от порядка 1

_ а m и n:

m-

ы— щ/ т—1 V

ТЖГ, , Бф) КsЩ2 + Ъ,s 1 + ... + ЪтУ2 + Ът ^ 5 ) алп2 + а, 5 "12 +... + ал -12 + а

и 1 п—1 п

Коэффициенты многочленов Д( -\|5) и Б( ) можно найти по частотной характеристике Ж (/Ю) аналитической передаточной функции (которая может быть получена экспериментально или на основе имеющихся временных характеристик). Вводя замену 5 —> /Ю и минимизируя

функционал ^^ Яе [ Б(/Ю) — Ж (/Ю) • Д(/Ю)]) , что отвечает методу наименьших квадратов

Ю

(МНК), определяем значения коэффициентов а и Ьу. Хотя частотная характеристика в общем случае — комплексная величина, тем не менее, хорошие результаты, в смысле регулярности оценок параметров, дает минимизация вещественной части функционала квадрата ошибки. В соответствии с частотной передаточной функцией, переписанной в виде:

-а0Ж(/ю) • (/ю) /2 = -ахЖ(/ю) • (/ю) /2 + -а2Ж(/ю) • (/ю)

п-2/

+ ...

+а „ -Ж(/ю) • (/ю) /2 + апЖ(/ю) - Ъ (/ю) /2 - ...- Ъп_х (/щ) /2 - Ът ,

при а = 1, составляется матрица регрессоров ^ размерностью [ NX (п + п +1)] , где N —

размерность вектора частот Ю , п — размерность вектора параметров а знаменателя, т+1 — размерность вектора параметров Ь числителя. Тогда, согласно классическому МНК, параметры передаточной функции вычисляются следующим образом:

[а !Ъ] = (Яе [Рт • Р ])-1 • Яе

где [а:Ъ] = [а .--ап Ъ0---Ъп] , Р = [р :р ]

Рт • -Ж(/ю) • (/ю ) -

Г п-1 п-2/

I (/ю ) /2 Ж (/ю ) (/ю ) /2 Ж (/ю )

п-1 п-2

(/ю ) ^ Ж (/ю ) (/ю ) /2 Ж (/ю )

С/ц) 12 ) жсМ ) (/ю ) ^ ЖС/ю ) Ж(/ю )

п-1 п-2 1

/) /2 ) (/Ю) /2 ) ... (jЮN) /2 ЖС^) )_

-(/ю )п

-( /ю1)

п -1

п п

-С/ю )п 2 -(/ю )

(МУ

-( /ю ) 1

— I

-1

п

) -(jЮN )

(jЮN ) -1]

Наиболее распространенным методом идентификации для указанного случая является известный метод Симою [2], позволяющий найти передаточную функцию целого порядка по снятой экспериментально кривой разгона. Осуществим модификацию этого метода для получения передаточной функции дробного порядка, что весьма актуально, применительно к теплотехнологиче-ским объектам.

Предположим, что исследуемый объект может быть описан линейным дифференциальным уравнением дробного порядка с постоянными коэффициентами:

Сп/2 х

а

с

п/ 2 + ап-1 -

С (п-1)/2 х

'А п-1)/ 2

+ ...+ а

С ^

1/ 2

+ а0 х = Ъ

С п/2Х

п/ 2

0

+ Ъп -1

С(п-1)/2 X

л

( п-1)/ 2

+ ... + Ъ

с1/2х

1 л

1/ 2

+ъох ,(1)

йг^' 2 1 Л1 2 0 п Сг

где а*,..., а*; Ъ° ,...,Ът — постоянные коэффициенты, х — отклонение регулируемой вели

чины, X — входное воздействие.

Приведем уравнение (1) к следующему виду:

Сп/2 х Лп-1)/2 х С1/2х Г . С п/ 2Х . Лп-1)/2Х . С1 2х 4

а„

Сг

п/ 2

а

п-1

Сг

-(п-ЦЙ

+ ... + а

Сг

-1/2-

х = К ( Ъ

п

I

С

п/ 2

Ъ

п-1

Сг ( п-1)/ 2

Ъ

Сг

1/ 2

- X

,(2)

где а ,.••, а„; Ъ,.., Ът — постоянные коэффициенты, К — коэффициент усиления, определяемый по формуле:

х

К = _ Уст

X

уст

Передаточная функция объекта Ж ( 5 ) = X ( 5 )|Л ( 5 ) , описываемого уравнением (2), быть представлена в следующем виде:

Ъ п5П/ 2

может

Ж (5) = КЖ (5) = К

+ъ 5п-1)/ 2 + ...+Ъ(51 2+1

п -1 1

ап5п/2 + а 5п-1)/ 2

п п-1

+...+а 51/2+1

2

2

2

2

2

2

2

2

2

2

2

2

2

п-1

*

*

В дальнейшем будем рассматривать нормированную передаточную функцию Ж (я) . Введем переменную V = , тогда передаточная функция объекта будет выглядеть следующим образом:

—, . Ь V т + Ь т-1 + ... + Ьь + 1

' 4 _т-1_1_

Ж (V) =

аьп + а„ V п 1 + ... + аь + 1

п п—1 1

Рассмотрим инверсную передаточную функцию объекта:

1 . . 1 а V п+ а .V п1 + ... + aV +1

__ _п_п—1_Л_

Ж (V) = = , .

" Ж (V ) Ь„Р т + Ьm-lV т-1 + ... + Ь? + 1

--1

Разложение Ж (V ) в ряд Тейлора в окрестности точки V = 0 имеет вид:

Ж4 (V) = 1 + C1V + Ср 2 + ... + ... , (3)

где коэффициенты разложения Ск называются согласно методу Симою площадями. При-

—-1

равнивая выражение для передаточной функции Ж (V ) и разложение (3), получим:

(1+ ая + а2я2 + ••• + апяп ) = (1+ Ья + Ь2я2 + ••• + Ьтят )(1 + Сх$ + С2я2 + ••• + Скяк + ). После раскрытия скобок и приведения подобных имеем:

к-1

Е

1 + а я + а я2 + ••• + а я = 1 + (Ь + С + (Ь + Ь С + С )я2 + ••• + (Ь+ Ь С )/ + С )/ + •

12 п 11 2 1 12 к ^ г к-г к

,2 , , ___, п. , г< \„ , /г. , ип . г< \„2 , , гг. , ип \„к , ^ \к

г=1

С целью определения неизвестных коэффициентов а1 ...ап; Ь1...Ьт составим линейную систему уравнений, приравнивая коэффициенты при одинаковых степенях 8 в последнем уравнении:

а = ь+С ;

а2 = Ь2 + Ь1С1 + С2 ;

а = Ь + Ь2СХ + ЬХС2 + С3; (4)

к-1

ак = Ьк + Ск + Е Ь'Ск-г .

г =1

Предположим, что идентифицируемый объект описывается передаточной функцией порядка п, а т = п-1, тогда, для получения значений коэффициентов а,...,ап; Ьх,..., Ьт необходимо составить систему из т + п уравнений.

Согласно классическому методу Симою идентификация математической модели объекта производится по кривой разгона, снятой экспериментально — Ых). Выбирая за начало отсчета точку Ы(0) и масштабируя кривую разгона путем деления на К, получим нормированную переходную

характеристику к (*) такую, что к (<^) = 1.

Введем вспомогательную функцию ф ) , определяемую формулой:

Л1 2

ф (*) = ^17ГI1 - к(*)) . (5)

Найдем изображение Ф ( я ) по Лапласу для оригинала ф ) :

Ф(я) = Ь {ф (*)} = лЛ -

1- Ж (я) 1 - Ж (я)

* '

Перейдем к переменной V :

ф(Щ)=1-Ш. (6)

V

Функцию Ф(и ) можно разложить в ряд Тейлора по степеням V в точке V = 0 :

Ф(и) = |0 + |и + |2и 2 + ... + |ки к + ... (7)

где коэффициенты разложения, называемые моментами вспомогательной функции, находятся по известной формуле:

1 С к Ф(и )

1 к к ! Си -

(8)

Моменты | могут быть определены непосредственно из функции ф (г ) , вычисленной по экспериментальной переходной характеристике. Вычисление операции полудифференцирования

С1/2

1/2 в выражении (5) будет рассмотрено ниже.

Изображение вспомогательной функции ф (г ) согласно преобразованию Лапласа есть функция:

Ф ( 5 ) = £ ф (г ) б?" 5г Сг

(9)

аргументом которой является комплексная переменная й, поэтому для вычисления моментов (8) производные Ф( к) (и ) =-(— целесообразно определить по правилу дифференцирова-

Си к

ния сложной функции:

СФ(5) ЭФ(и) Си

С5 Эи С5

Тогда рекуррентная формула для нахождения к-й производной будет иметь вид:

Ф(к) (и) = ЭФ(к-1)(и) = СФ(к-1) (.у»и _ Эи ё5

Учитывая, что и = и С^Си = 2*//5 имеем:

ЭФ(и) =2 ^ Ф(5) .

Эи ё5 '

Э 2Ф(и ) = 2 С Ф(5) + ^ С2Ф(5) .

Эи 2 С5 С52

^=12 /15^+85 2 С^

Эи 3 С5 С5

ЭкФ(и) Ф(к-1) (5)

Эи к С5 '

где производные к-го порядка функции (9) имеют вид:

С к Ф ( 5 )

к

к ш I V I

ш' 4' = | ( -г)кф ( г) Сг . (10)

С5

Искомые моменты (8) найдем, подставив в (10) 5 = 0 , с учетом свойства преобразования Лапласа о дифференцируемости функции-оригинала:

и =0

0

I = / Ф (г )Сг;

0

11 = 2 ^ (-г )Ф (г С;

1/ 2

1 Г ~ С - ^

I = —2 (-г )ф (г )Сг + 4 — (-г )2 ф (г )Сг ; (11)

2 2! [ I Сг I у

1 Г С17 2 - , ,2 , ^ 0 С 3 7 2 - ^

I = -12- (-г )ф (г )Сг+ 8- (-г )3ф (г )Сг .

3 3! [ Сг1 2 | Сг 3 / 2 |

Установим теперь связь между коэффициентами | и С, для чего преобразуем формулу (6):

и • Ф (и ) = 1 - Ж (и )

или

Ж (и ) = 1 -и • Ф (и ) .

Откуда следует,

(1-и • Ф (и )) • Ж 4 (и ) = 1 . В последнее выражение подставим разложения (3) и (7):

[1 - |0и - ци 2 - | 2и 3 ... - к - ...] (1 + Си + Си 2 + ... + Си к + ...) = 1 .

Л.0и 1 2и ... 1к -1и . . 1 С1и 1 С2и 1 ... 1 Ски

Раскрывая скобки и приводя подобные члены, получим:

Г к -2

1 + ( С1 - I )и + ( С2 I С1 - 11 )и 2 + . " + Ск - I -1 - Е I Ск-1

- к Г~к -1 ¿^Ъ к -1-1

V 1=0 у

и 2 + ... = 1

Приравнивая коэффициенты при одинаковых степенях и , составим линейную систему уравнений:

С - I = 0;

С2 - |0С1 - 1 = 0;

к -2

Ск - 1к -1 - Е 1Ск-1--=0

1 =0

решая которую, найдем рекуррентные формулы для вычисления коэффициентов С:

С1 = 1 ; С2 = |0С1 + 1;

Ск=I -1+Е I с

к -2

к -1-1 .

Далее, подставляя значения коэффициентов Ск в систему (4), определим искомые параметры передаточной функции а,. • •, ап; Ъ,-., Ът •

Вычисление производных дробного порядка, кратного ^ , переменных в выражениях (5), (11) может быть основано на определении дробной производной Летникова [3]. Дробная производная порядка V от непрерывной функции является также непрерывной функцией следующего вида:

I=0

— 1'(* В' f ( *) =

—т~-Г ' — Г(* -Т ' 1 (т ) —Т , -^ <V < +1

Г (1 -V ) — 0

Для частного случая V зависимость примет вид:

или ( )

* 1

Г —

В2 / ( * ) = 1 • — /

( )

( )

* - т

1 1 1 2 /'^_П7->2 _ПГ2

Б2 / * = В Б 21 * = В I 2 / * , где 1 — операция полуинтегрирования.

Операция дифференцирования В целого порядка выполняется по одному из известных разностных алгоритмов. Операция полуинтегрирования осуществляется на основе численного интегрирования

121 (* ) = . Г/(Т ) ^

Ф 0 '

V - т

После перехода к решетчатой функции, в результате чего в случае применения метода прямоугольников на к-м шаге, получим

-2 Г П т к-1 / [Т ] 2 [ ] -г- е -

I / кТ =

где Т0 — шаг по времени. П г =0 (к - г ) Т0

Проиллюстрируем возможности предложенной модификации метода Симою для идентификации объекта передаточной функцией дробного порядка на конкретном примере. Рассмотрим задачу торцевого нагрева полубесконечного тела, описываемого одномерным односвязным уравнением второго порядка параболического типа:

дТ (г, *) Э2Т (г, *)

-= 0.05-.

д* д2 г

Начальная температура тела 0° С , температура на торце полубесконечного тела — 100° С . С помощью метода сеток построим кривую нагрева тела на глубине Г = 1 м (рис. 1).

10 20 30 40 50 50 70

Рис. 1. Кривая нагрева сечения стержня для Г = 1

Рис. 2. Результаты аппроксимации

Сравним результаты идентификации объекта, проведенные классическим методом Симою с получением передаточной функции целого порядка и модифицированным методом, дающим описание в форме передаточной функции дробного порядка.

В результате идентификации объекта передаточной функцией дробного порядка модифицированным методом Симою, изложенным выше, имеем:

Ж (я) = 100

902.84я + 7.53я05 + 1

1857.74^ 5 + 906.94я + 10.06я0 5 + 1

Применяя к этому же примеру метод идентификации передаточной функцией дробного порядка по ВЧХ с помощью МНК, имеем:

Ж (я) = 100

0.45я -1.09/5 + 0.97

3.69 -10-5 я15 + 8.9 -10-3 я + 0.12/5 + 1

При идентификации объекта передаточной функцией целого порядка, используя средства модуля ident пакета МаЙаЪ, получаем:

м

TXW ч n ™ U.UUU25U9s + U.UU1U88s + U.UU244s + U.UU2485

W3 (s) = 1UU-----

3 s3 + U.8737s2 + U. 1447s + U.UU27U6

На рис. 2. сплошной линией показаны исходные данные, полученные путем решения крае- вой задачи теплопроводности для рассматриваемого примера методом сеток, штриховой линией показан результат аппроксимации модифицированным методом Симою (СКО составляет 1.99), штрихпунктирной — показан результат аппроксимации по ВЧХ (СКО составляет 6.26), пунктирной

— показан результат аппроксимации передаточной функцией целого порядка (СКО составляет

18.37). Как видно на рис. 2., при использовании модели целого порядка наблюдается существенное

отклонение ее переходной характеристики от исходной кривой при значениях времени, близких к

установившемуся режиму. Обе модели дробного порядка (рис. 2, штриховая и штрихпунктирная

линии) довольно точно описывают характер переходного процесса исследуемого объекта на протяжении всей длительности эксперимента.

Необходимо отметить преимущество изложенного выше метода идентификации в услов и-

ях реального производства, обусловленное сравнительной простотой построения кривой разгона, в

отличие от идентификации по ВЧХ, так как построение частотной характеристики реального объекта сопряжено со значительными временными и энергетическими затратами, в особенности для

весьма инерционных теплотехнологических объектов автоматизации.

Работа выполнена при финансовой поддержке ФЦП «Научные и научно-педагогические

кадры инновационной России» на 2009-2013 годы. Госконтракт № 14.740.11.0591 от 05

октября

2010 г.

Список литературы

1. Мишунин, В.В. Идентификация передаточных функций электрических печей цилиндрической формы // Передовые технологии в промышленности и строительстве на пороге XXI века: Сб. трудов. - Белго- род, 1998. - Ч. 1. - С. 1025-1031.

2. Симою, М.П. Определение коэффициентов передаточных функций линеаризованных звеньев систем регулирования// Автоматика и телемеханика, 1957. — № 6, — с. 514-527.

3. Бабенко, Ю.И. Метод дробного дифференцирования в прикладных задачах теории тепломассо- обмена. — СПб.: НПО «Профессионал», 2009. — 584 с.

MODELING OF TECHNOLOGICAL HEATING OBJECTS IN THE CLASS OF FRACTIONAL IRRATIONAL TRANSFER FUNCTIONS

E.B. KARIKOV V.V. MISHUNIN2 V.G. RUBANOV Y.A. GOLTSOV

1) Belgorod Shukhov State Technology University

2) Belgorod National Research University

e-mail: rubanov@intbel.ru

This article describes methods for identification of distributed pa- rameter systems using frequency and time characteristics. For identifi- cation by the time characteristics was used a modified Simoyu identifi- cation method. Shows the results of the identification of different me- thods on the example of solving the problem of heating a semi-infinite body.

Keywords: identification, fractional irrational transfer function, the method of least squares, Simoyu.