http://habrahabr.ru/post/218297/

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

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

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

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

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

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

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);

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

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

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

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

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

omega = (2*%pi*%i*frequency);
tmp = sp.*omega;
out = real(fft(tmp,1));
speed_temp = out(1:N/2);

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

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);

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

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 );

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

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 ;

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

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;

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

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

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--");

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


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