Идентификация быстрых термических процессов

Недавно мне удалось завершить часть работы по очень интересному проекту в ФТИ им. Иоффе и получить достаточное количество экспериментальных данных, для того чтобы поделиться с Вами.
Физики из СПб ФТИ им. Иоффе занимаются выращиванием нитрид галлиевых полупроводниковых структур, которые обладают неплохими показателями скорости носителей заряда при переходе и большим коэффициентом теплопроводности. Процесс роста такой структуры проходит при температуре 1000 С (1273 К) и атмосферном давлении. Все происходит в специальной камере, находящейся в герметичной зоне. При выращивании структуры весь объем реактора и герметичной зоны заполняется азотом. В процессе роста структуры подложкодержатель вращается с частотой один раз в секунду. Такие операции относятся к быстрым термическим процессам, скорость изменения температуры в которых варьируется от нескольких единиц до сотен градусов в секунду.



Моей задачей было управление температурой графитового подложкодержателя при помощи индуктивного нагрева.
Технические характеристики установки выглядят следующим образом. Для измерения температуры используется лазерный пирометр, снимающий данные в центре графита. Частота съема информации 10 раз в секунду, шаг измерения 1 градус. Значение мощности передаваемой графиту полагается прямо пропорциональным мощности на индукторе. У генератора управляющего индуктором имеется цифровой выход, по которому передается напряжение, ток и мощность.

Для начала нужно было настроить регулятор температуры, чтобы в процессе роста не было сильных колебаний. С этой задачей справились достаточно быстро, но наше решение давало качественный результат только на высоких температурах. Для других процессов требовалось менять коэффициенты.
Так как эта работа вписалась в тематику моей кандидатской, то захотелось написать какой-нибудь хитрый алгоритм управления на основе анализа модели термических процессов. Для начала ознакомился с законом Стефана — Больцмана, мощность излучения абсолютно чёрного тела прямо пропорциональна площади поверхности и четвёртой степени температуры. С учетом конвекции можно записать уравнение для термических процессов в точке съема температуры

где T — температура, P — мощность, Tc — температура окружающих объектов, которые нагревают графит собственным излучением, Ta — температура газа возле точки съема информации, который осуществляет конвекцию, B, A1, A2 — коэффициенты, которые необходимо идентифицировать. Для упрощения заменим все составляющие уравнения, которые мы можем считать константными в установившемся процессе и будем рассматривать следующее уравнение

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


Если мы получим оценку скорости температуры, то можно будет составить регрессионную модель с известными температурой и мощностью управляющего индуктора, и посчитать коэффициенты в уравнении.
Я уже достаточно давно использую для работы Scilab и в этот раз решил не изменять себе.Для этой задачи я выбрал взятие производной в образе Фурье. Но для начала работы с реальными данными, нужно интерполировать измерения для получения равномерной оси времени.

<code class="matlab">xx = linspace(0, round(time($)), round(time($))*fs+1)'; //New time axes' d=splin(time,temp,"monotone"); [int_temp,int_temp_diff] = interp(xx, time, temp, d); </code>

Стоит отметить, что переменная «int_temp_diff» будет содержать в себе данные скорости, но выглядят они весьма не лицеприятно.


Так что перейдем к работе с данными в образе Фурье. Нам нужно будет нарастить дополнительный хвост данных, чтобы на конце не было разрыва, зеркально отобразив график данных температуры.


Потом задаем ось частот и делаем дискретное преобразование Фурье.

<code class="matlab">fs = 10;           //Sample frequency in = [int_temp;int_temp($:-1:1)]; N = length ( in ); frequency = [0 : (N-1)] / N * fs; frequency (N/2+1 : $) = frequency (N/2+1 : $) -fs; frequency = frequency';  sp = fft(in);    //Fast Fourier Transform </code>

На частоте 1 герц виден явный пик, это связано с тем, что подложкодержатель вращается с этой частотой. Другие частоты пробились из-за того что графит нагрет неравномерно, и при вращении температура может увеличиться и уменьшиться несколько раз.
Чтобы взять производную в образе Фурье нужно просто домножить на. После этого делаем обратное преобразование Фурье.

<code class="matlab">omega = (2*%pi*%i*frequency); tmp = sp.*omega; out = real(fft(tmp,1)); speed_temp = out(1:N/2); </code>

Это конечно немного лучше, чем то что возвращает функция взятие производной через кубические сплайны, но идентификация с этим выдает некачественные результаты. Придется прибегнуть к фильтрации сигнала.
По совету моего друга Юры, который занимается настройкой корабельных систем управления, было решено использовать окно Блэкмана-Харриса для вырезания лишних частот. Так как в процессе роста графит вращается с частотой 1 раз в секунду, все колебания с частотой выше этой нам не интересны. Поэтому стоит вырезать все колебания выше 1 [Гц]:

<code class="matlab">cf = 1; //Cutoff frequency k = round(cf*N/fs); L = 2*k+1;  a0=0.35875; a1=0.48829; a2=0.14128; a3=0.01168;  for j=0:L-1,     w(j+1) = a0 - a1*cos(2*%pi*j/(L-1))+a2*cos(4*%pi*j/(L-1))-a3*cos(6*%pi*j/(L-1)); end h = zeros(frequency); for j=1:1:k+1     h(j) = w(k+j); end   h([$:-1:$-k]) = h([1:1:k+1]); omega = (2*%pi*%i*frequency); omega = omega.*h; tmp = sp.*omega; out = real(fft(tmp,1)); speed_temp = out(1:N/2); </code>

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

<code class="matlab">constant (1:length(int_power)) = 1; norm_4_temp = norm( int_temp.^4 ); norm_int_temp = norm( int_temp ); norm_int_power = norm( int_power ); norm_speed_temp = norm( speed_temp ); norm_constant = norm( constant ); </code>

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

<code class="matlab">LSM = [ int_temp.^4/norm_4_temp , int_temp/norm_int_temp , int_power/norm_int_power , constant/norm_constant ] ; a0 = pinv(LSM) * speed_temp / norm_speed_temp ; </code>

И после этого необходимо вернуться к оригинальным значениям вектора коэффициентов.

<code class="matlab">a0(1) = a0(1)*norm_speed_temp/norm_4_temp; a0(2) = a0(2)*norm_speed_temp/norm_int_temp; a0(3) = a0(3)*norm_speed_temp/norm_int_power; a0(4) = a0(4)*norm_speed_temp/norm_constant; </code>

В результате получились следующие значения a0 = [ — 1.073D-12, — 0.0029096, 0.0004617, 2.0969723 ], что достаточно близко к результату наших зарубежных коллег, хотя у них используется ламповый нагрев.
Вот теперь можно приступать к проверке результатов на полученной модели. Воспользуемся для этого пакетом виртуального моделирования Xcos Scilab. Так как у нас есть дифференциальное уравнение, описывающее термический процесс в реакторе, и все коэффициенты к нему, то соберем следующую схему.


Блоком «From Workspace» подаем экспериментальные данные мощности, предварительно для этого блока создав структуру V. Задаем частоту сэмплирования и время начало моделирования в блоке «Clock». Длину выходной структуры данных задаем равной входной. Конечное время моделирования должно равняться времени эксперимента. Составленную модель запускаем в тексте программы и сравниваем графики эксперимента и моделирования.

<code class="matlab">V.time = xx; V.values = int_power; importXcosDiagram("D:\PhD\Term_model.zcos"); xcos_simulate(scs_m,4); plot(time,temp); plot(A.time,A.values,"r--"); </code>

Вот еще тоже самое, но для гармонического возмущения с изменяющейся частотой.




Огромное спасибо научному руководителю Arastas за неоценимую помощь и участие в подготовке материала. Благодарю Заварина Евгения за предоставление доступа к установке и программную реализацию всех экспериментов.

Источник: habrahabr.ru/post/218297/