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

3

ЭНЕРГОМОНИТОРИНГ И ЭНЕРГОСБЕРЕЖЕНИЕ

МЕТОДИЧЕСКИЕ ПОГРЕШНОСТИ ПАРАМЕТРИЧЕСКОЙ ИДЕНТИФИКАЦИИ МОДЕЛЕЙ ТЕПЛОПЕРЕНОСА В НЕСТАЦИОНАРНОЙ ТЕПЛОМЕТРИИ Н.В. Пилипенко

Предложен универсальный для различных преобразователей теплового потока (ПТП) метод оценок точности результатов решения, а также оптимального планирования реализации граничной обратной задачи теплопроводности (ОЗТ) по восстановлению плотности входящего в ПТП теплового потока д(т). Метод заключается в построении совместных доверительных областей (СДО) или интервалов (СДИ) для оптимальных оценок параметров, реализующих кусочно-линейную сплайн-аппроксимацию д(т). Для построения СДО или СДИ использована их известная связь с ковариационной матрицей (информационной матрицей Фишера) ошибок указанных оценок. Элементами этих матриц являются комбинации функций чувствительности измеряемых температур ПТП к искомым параметрам д(т). Приведены результаты комплексных исследований особенностей реализации предложенной методологии нестационарной тепло-метрии и оценок ее погрешностей применительно к некоторым распространенным ПТП.

Восстановление нестационарного теплового потока д(т) с помощью различных ПТП связано, как правило, с рядом методических погрешностей, основные из которых представлены на рис.1. При этом можно выделить общие погрешности, присущие как нестационарной, так и стационарной теплометрии, вызываемые неоднородностью измерительной среды - искажающим влиянием ПТП на теплообмен и температурное поле объекта исследований. Эти источники погрешностей детально исследованы в работах [1-4] и в данной работе рассматриваться не будут.

Принципиальной особенностью нестационарной теплометрии является то обстоятельство, что в общем случае приходится преодолевать тепловую инерционность ПТП расчетным путем - восстанавливать значение входящего в ПТП теплового потока д(т) по измерениям температур или их перепадов в ПТП. Возникающие при этом погрешности, вызываемые процедурой восстановления д(т), могут быть отнесены к основным погрешностям нестационарной теплометрии. Они являются весьма существенными, а в ряде случаев - определяющими. Сложности решения проблемы устранения и оценивания основных погрешностей определяются следующими обстоятельствами:

• задача восстановления д(т) относится к некорректно поставленным обратным задачам теплопроводности (ОЗТ), и ее решение вносит существенный вклад в погрешность нестационарной теплометрии;

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

Для упрощения решения задачи нами был принят и реализован ряд требований и допущений к конструкции ПТП, математическими моделями теплопереноса в них, к характеру изменения восстанавливаемого теплового потока д(т) и к методам его восстановления, а именно:

• требование обязательного измерения либо температуры рабочей поверхности ПТП, либо перепада температур между нею и другой точкой ПТП. Его выполнение переводит граничную ОЗТ в разряд псевдообратных, что снижает сложность алгоритмов

восстановления д(т) и общий уровень основных погрешностей теплометрии, а также естественным образом предоставляет информацию о моменте начала действия д(т), принципиально необходимую для решения ОЗТ; • допущение, что априорные сведения о характере измеряемого теплового потока д(т) позволяют выполнить его кусочно-линейную В-сплайн-аппроксимацию и на

каждом ее г-ом участке выделить вектор искомых параметров ( = 1,2,3,...), представляющий д(т).

Погрешность решения ПЗТ

Погрешности неадекватности

модели теплопереноса

Структурные

Системати- Случайные

ческие

Параметрические

Алгоритма минимизации

Рис.1. Структура методической погрешности нестационарной теплометрии

на основе одномерных ПТП

Это дает возможность использовать для восстановления д(т) метод параметрической идентификации ПТП, в частности, алгоритм фильтра Калмана (ФК) [5], а также разделить основную погрешность нестационарной теплометрии на две составляющие -погрешность аппроксимации д(т) и погрешность параметрической идентификации.

Предположим, что, исходя из анализа информации о физической природе и термогазодинамических свойствах среды, с которой находится в теплообмене рабочая поверхность ПТП, априорная сплайн-аппроксимация д(т) может быть выполнена с необходимой точностью. Тогда в соответствии с рис. 1 основной погрешностью будет погрешность параметрической идентификации ПТП - процедуры последовательного определения оптимальных оценок Qz е вектора искомых параметров Qz на 2 = 1,2,3, к участках сплайн-аппроксимации д(т). Вектор Qz может быть представлен в виде

(? г =

да, дЬг

где да и дЬ - значения потока в начале и в конце участков сплайн-

аппроксимации.

Для нахождения оптимальных оценок Qf на г-ом участке нами используется алгоритм фильтра Калмана, а в качестве исходной информации - 1 измерений температур в ПТП на этом участке, составляющих (т х 1)-вектор измерения Ук, где к = 1,2,3,..., 1 - дискретное время.

Фильтр Калмана в силу особенностей своей практической реализации, по сути, является рекуррентной процедурой обобщенного метода наименьших квадратов (МНК), минимизирующей по вектору искомых параметров Q функцию невязки ф(

Т

Ф(( )=!

к=1

Ук - ^

(к) [4 - Ук (к)

(1)

где Ук - вектор измерений температуры ПТП, включающий вектор ек случайных погрешностей измерений; Ук (1к) - модельные (расчетные) значения вектора измерений; к - дискретное время; Т- знак транспонирования.

Для линейных МНК, когда вектор Ук и его модельные значения Ук () линейно зависят от вектора искомых параметров Q, известны формулы [6] для оптимальных

оценок Qf и ковариационной (2 х 2) матрицы ее ошибок Р1 Qf 1:

к=1

дQ

Qt = Р^М _

^ <21

Р^) = Ц НТК-1 нк л 1

ук = Р Q1 нТ • у ;

к=1

к=1

(2) (3)

где Нк - (т х 2) -матрица функций чувствительности, которая имеет вид

дУ1к дУ1к

Нк =

ГдУк >

дQ

дда ддь

дУтк дУтк

дда ддь

и 11к и 12к

^ т1к ^ т 2к

Ql

Ql

. (4)

Матрицу Нк составляют функции чувствительности и ук, которые количественно характеризуют чувствительность /-ого измерения у к в к-й момент времени к изменению у-го искомого параметра да (( = 1) или дь (7 = 2). По определению,

ипк =

дУк

дда

иг 2к =

дУк

Ql

ддь

(5)

Ql

В теории МНК показано, что ковариационная матрица Р^ ) является характеристикой точности оценок Q1. В рассматриваемом случае она имеет вид

Р11 Р12

Р1 Qí ) =

(6)

к

21

22

Ее диагональные элементы р11 и р22 являются дисперсиями оценок даг и цъ г

соответственно, а элементы р12 = р21 характеризуют взаимную корреляцию этих оценок.

В формуле (3) разделим влияния уровня а2 шума в измерениях и других значимых особенностей (факторов) нестационарной теплометрии, которые отражены в матрице функций чувствительности Ик, и получим:

РШе} = а2-IIИТ - И

V1

Т

к ±± к

к=1 У

(7)

Выражение в скобках является (2 х 2) -матрицей Грама Ле для системы векторов функций чувствительности (она же - аналог информационной матрицы Фишера). В рассматриваемом случае

ле = 1 ИТ - Ик =

к=1

IIЦип1к - и

к=1 I =1 к=1 1=1

1 п 1 п

II Рп 2 к - ип1к II и

п2к

2

п2к

(8)

к=1 1=1 к=1 1=1 Введем понятие характеристической ковариационной матрицы:

р = Л -1 = Р11 Р12 Г1 - - Р — •

Р 21 —22

Тогда выражение (7) для ( } примет вид

( ) = а2 - р . Из (8) и (9) следует соотношение:

р(( У-1 р= Л£,

К У а а

(9)

(10)

которое будет использовано в дальнейшем.

В качестве иллюстрации приведем матрицы функций чувствительности Ик и Грама для двух вариантов измерения температур в ПТП. При измерениях температуры в одной точке ПТП на поверхности (п=1) (вектор измерений - Ук = t1k +ек) матрицы Ик и Ле имеют вид

Ик = Р

1к Р 2к

Ле =

(1х2)

IV 12к IV 1к - Р

к=1 к=1

IV 2к - Р1к Р 2к

(11)

к=1 к=1

При двух измерениях температуры в ПТП (п=2), в частности, температуры поверхности t1k и перепада температур, например, Д^1-6 )к = t1k -16к , (2 х 2)-вектор измерения имеет вид

П =

У1к +е1к tlk

У 2к + е2к ^ (1-6 ),к

+ ек,

а матрицы функций чувствительности Ик и Грама Ле

Ик =

Р 11к Р 12к Р 21к и 22к

(12)

(13)

е

I

л =

+ и221к) К • и12к + и

21к ' и22 к )

к=1 к=1 I ((12к • и 11к + и22к • и21к ) I ((121к + и22к )

(14)

к=1 к=1

Нами проанализированы все источники погрешности параметрической идентификации д(т). Однако в данной работе основное внимание уделено двум всегда имеющим место погрешностям - влиянию случайных погрешностей (шума) в измерениях температур ПТП и погрешности, вызываемой топологией функции невязки Ф(й).

Достоинство предложенной методологии нестационарной теплометрии и оценок ее основных погрешностей заключается в том, что указанные составляющие погрешности количественно отражены в ковариационной матрице р(й ), и в то же время эти две

погрешности из формулы (3) для р(й) можно четко разделить на два сомножителя -

дисперсию а2, определяющую уровень шума в измерениях температур ПТП, и характеристическую ковариационную матрицу р, которая через матрицу функций чувствительности Нк выражает основные особенности ПТП как средства измерения. Остановимся на этом более детально.

Ковариационная (2 х 2)-матрица р(й) ошибок оценок й в соответствии с (3) получается путем обращения (2 х 2)-матрицы Грама, элементами которой являются функции чувствительности измеряемых температур в ПТП Ук к составляющим ца и дь искомого вектора й на 2-ом участке его аппроксимации.

Ковариационная матрица й ^ в соответствии с (7) через дисперсию а2 погрешностей в измерениях учитывает случайные погрешности средств измерения температур и соответствующей регистрирующей аппаратуры, а через характеристическую ковариационную матрицу р - все другие значимые факторы процесса нестационарной

теплометрии, такие как вид, количество и характер размещения в ПТП непосредственно измеряемых температур или их перепадов; характер переходного процесса в ПТП,

включая начальное распределение температур Т0; закон изменения восстанавливаемого д(т); участок г наблюдаемого переходного процесса; количество измерений 1 вектора измерений Ук на этом участке (к = 1, 2, ... 1).

Как было отмечено ранее, основной погрешностью в нестационарной тепломет-рии после учета остальных, представленных на рис. 1, является погрешность параметрической идентификации. Для ее оценки нами использовался подход, основанный на построении совместных доверительных областей (СДО) и интервалов (СДИ) оценок искомых параметров [5]. Рассматривалось два случая: единичного (скалярного) ш=1 и векторного ш=2,3,... измерений температур или их перепадов в ПТП. Для случая скалярного (ш=1) измерения Ук = ук в соответствии с [6] уравнение СДО для оценок ^ и двух параметров имеет следующий вид:

й - 4 )г • л {е - 4 )=^2 • 2 • 1 - 2), (15)

где Fv (2,1 - 2) - табличное значение квантиля распределения Фишера для случая 1 измерений на каждом г-ом участке; V - величина доверительной вероятности; £2 -

выборочная дисперсия шума в измерениях, определяемая непосредственно, исходя из экспериментальных данных, по формуле:

о

1

52 =

1 - 2

•I

к=1

ук - ш

Учитывая ограниченное количество 1 измерений Ук на участках аппроксимации,

для оценочных расчетов в качестве 52 может быть взята дисперсия а2, соответствующая паспортным данным о средствах измерения температуры ПТП.

Для случая векторного измерения (т = 2, 3,...) в соответствии с [6] СДО для оценок тех же двух параметров да и Чь имеет следующий вид:

б - &] • Л • [ б - &] =

2

1 - 2

Ф • ^ (2 1 - 2)

шт VV ' !'

где Фшт = Ф\бА - минимальное значение функции невязки при б = .

(16)

Доверительный эллипсоид (16) отличается от (15) операцией вычисления минимального значения функции невязки Фшт взамен выборочной дисперсии 52 по следующей формуле, вытекающей из (7):

Ф ■ = а21

шт /

к=1

Ук - Ук [б.

У - У

кк

(17)

Для рассматриваемого случая СДО могут быть построены методом канонического анализа эллипсоидов (15) и (16). Очевидно, что они в двухмерном случае имеют форму эллипсов в пространстве двух искомых параметров.

Пусть получены оценки =

Чаи

Чье

локального вектора искомых параметров б на

2-ом участке кусочной сплайн-аппроксимации ч(т), и на их основании по формулам (8) или (11), (14) рассчитана матрица Грама Ае, которую запишем в виде:

Ае =

а11 а12

а21 а 22

Тогда уравнения (4.17) и (4.18) совместной доверительной области принимают

вид:

где

[(Ча - Фае )(Чь - €)]

а,

а

а12 Ча - ФсЛ

а 22 _ _Чь - Фы _

В =

252(2,1 - 2)

2 Фшп Ру(2,1 - 2)

1 - 2

= В.

■ для скалярного измерения (ш = 1);

■ для векторного измерения (ш > 1).

(18)

Проведя преобразования, получим следующее уравнение кривой второго порядка (эллипса) в пространстве двух искомых параметров Ча и Чь :

а11 •(Ча - фа1 )2 + 2а12 (Ча - €а1 ) • (Чь - ^Ы )+ а22 • (Чь - ^Ы - В = 0 (19)

Все коэффициенты уравнения (19), в том числе и свободный член В, либо известны, либо определяются, исходя из условий проведения эксперимента, по приведенным выше зависимостям. Для рассмотрения особенностей полученной области удобным является способ построения эллипса по точкам, заключающийся в следующем: задаваясь значениями одного из параметров, например, Чь , можно получать из (19) квадратные

Т

уравнения относительно отклонения второго параметра, например, (да - ), которые в системе координат (да - ) и (дь - ) дают по два значения на кривой эллипса. Решения этих уравнений имеют вид

( £ ) - а12 - §ы )± V а122 - (Чь - §ы )2 - а11 • а22 -(^ь - §ы )2 -

Ща - ^Ра« /1,2 ---(20)

или

( д 2 - - а12 • (а - С 2±У а122 • (а - С )2 - а22 • [а11 • (а - С )2 - В

Щь - ? /1,2 --

а22

(21)

Таким образом, для случая любой размерности т вектора измерений Ук и размерности г =2 локального (для г-го участка сплайн-аппроксимации) вектора Q могут быть

достаточно просто построены двухмерные СДО оценок 0«, имеющие наглядную форму эллипсов в пространстве да и дь.

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

^с (о - ^с (г), (22)

где х1 (г) - квантиль х2 -распределения для г параметров (в рассматриваемом случае г=2); Q0 - задаваемые априорно «истинные» значения вектора параметров.

В соответствии с формулой (10) уравнение (22) может быть преобразовано к виду

'- 0?0Т А{?- 0?о 1- ° х1 (г). (23)

Уравнение (23) отличается от (15) и (16) значениями правых частей. Очевидно, что с учетом этих отличий для построения СДО или СДИ при оптимальном планировании могут быть использованы все предыдущие зависимости. Однако предварительно должны быть выполнены следующие этапы планирования.

1. Исходя из априорной информации об объекте и условиях проведения тепло-метрических исследований:

1.1. выбрать тип и параметры ПТП;

1.2. путем сравнительного анализа динамических характеристик ПТП и теплообмена объекта со средой убедиться в необходимости решения граничной ОЗТ по восстановлению д(т);

1.3. на основе априорных представлений о д(т) обосновать принципиальную возможность его кусочно-линейной сплайн-аппроксимации и выбрать вектор искомых параметров ?.

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

убедиться в возможности получения оценок ? «, близких априорно задаваемым «истинным» значениям ? - .

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

3. Задаемся требуемыми характеристиками точности идентификации - значениями СДИ для каждой составляющей б или формой и размерами их СДО (в случае г < 2).

4. Задаемся априорными истинными значениями б0 вектора идентифицируемых параметров б .

5. Задаемся комбинацией ранее перечисленных значимых факторов эксперимента и их числовыми выражениями. По ДРМ ПТП с учетом указанных факторов и величины

б0 вычисляем функции чувствительности и]1к (/=1,...,да;у=1,...,г; к=1,...,Ы), а далее по формулам (4) или (13) - матрицы чувствительности Йк . В указанных формулах вместо

б подставляем б.

6. В соответствии с формулой (8) на основе полученных матриц Йк и количества измерений 1 на рассматриваемом 2-ом участке получаем матрицу Грама.

7. Для случая г < 2 можно ограничиться полученной матрицей Грама и на ее основе построить СДО и СДИ для оценок б, аналогичные описываемым формулами

(20). Однако в них, в соответствии с видом доверительного эллипсоида (23), необходимо подставлять значение

В = а2х1 (2), (24)

где а2 - априорная оценка дисперсии помехи в измерениях температур ПТП.

8. Путем сравнения полученных значений СДО с заданными (п. 3) произвести вывод о правильности выбора значимых факторов эксперимента. При отрицательном результате процедуру планирования необходимо повторить до удовлетворения заданных требований.

В качестве примера планирования эксперимента приведем построение СДО простейшего ПТП в виде вспомогательной стенки при измерении температуры ее поверхности ^(т) со среднеквадратичным отклонением а = 1 °С; при условии

б =

Ча 105

Чь 6-105

Вт/м2.

Размеры и свойства материала ПТП примем следующие: толщина к = 5 -10 м;

Вт м2

теплопроводность X = 15-; температуропроводность а = 3,8 -10-6 —. Пользуясь

м • К с

приведенной выше последовательностью, получим, что погрешности на первом участке

3 2 3 2

сплайн-аппроксимации равны АЧа = 1 -10 Вт/м , Ачь = 3.5 • 10 Вт/м , а совместная

доверительная область имеет вид рис. 2. Построив СДО для всех участков сплайн-аппроксимации, получим суммарную погрешность восстановленного искомого потока

Ч(т).

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

и

Рис. 2. Совместные доверительные области параметров qa и qb

Литература

1. Ярышев Н.А. Теоретические основы измерения нестационарных температур. Л.: Энергоатомиздат, Ленинградское отделение, 1990. 255 с.

2. Геращенко О.А. Основы теплометрии. Киев: Наукова думка, 1971. 191 с.

3. Точность контактных методов измерения температуры / Под ред. Гордова и др. М.: Изд. стандартов, 1976. 232 с.

4. Коротков П.А., Лондон Г.Е. Динамические контактные измерения тепловых величин. Л.:Машиностроение, 1974. 224 с.

5. Симбирский Д.Ф. Температурная диагностика двигателей. Киев: Техника, 1974. 208 с.

6. Химельблау Д. Анализ процессов статистическими методами. М.: Мир, 1973. 957 с.