Недавно мне удалось завершить часть работы по очень интересному проекту в ФТИ им. Иоффе и получить достаточное количество экспериментальных данных, для того чтобы поделиться с Вами.
Физики из СПб ФТИ им. Иоффе занимаются выращиванием нитрид галлиевых полупроводниковых структур, которые обладают неплохими показателями скорости носителей заряда при переходе и большим коэффициентом теплопроводности. Процесс роста такой структуры проходит при температуре 1000 С (1273 К) и атмосферном давлении. Все происходит в специальной камере, находящейся в герметичной зоне. При выращивании структуры весь объем реактора и герметичной зоны заполняется азотом. В процессе роста структуры подложкодержатель вращается с частотой один раз в секунду. Такие операции относятся к быстрым термическим процессам, скорость изменения температуры в которых варьируется от нескольких единиц до сотен градусов в секунду.
![](https://habrastorage.org/r/w780q1/getpro/habr/post_images/1b0/9cc/cb4/1b09cccb471fa6f01aa7430f3d6f428b.jpg)
Моей задачей было управление температурой графитового подложкодержателя при помощи индуктивного нагрева.
Технические характеристики установки выглядят следующим образом. Для измерения температуры используется лазерный пирометр, снимающий данные в центре графита. Частота съема информации 10 раз в секунду, шаг измерения 1 градус. Значение мощности передаваемой графиту полагается прямо пропорциональным мощности на индукторе. У генератора управляющего индуктором имеется цифровой выход, по которому передаются значения напряжения, тока и мощности.
Для начала нужно было настроить регулятор температуры, чтобы в процессе роста не было сильных колебаний. С этой задачей справились достаточно быстро, но наше решение давало качественный результат только на высоких температурах. Для других процессов требовалось менять коэффициенты.
Так как эта работа вписалась в тематику моей кандидатской, то захотелось написать какой-нибудь хитрый алгоритм управления на основе анализа модели термических процессов. Для начала ознакомился с законом Стефана — Больцмана, мощность излучения абсолютно чёрного тела прямо пропорциональна площади поверхности и четвёртой степени температуры. С учетом конвекции можно записать уравнение для термических процессов в точке съема температуры
где T — температура, P — мощность, Tc — температура окружающих объектов, которые нагревают графит собственным излучением, Ta — температура газа возле точки съема информации, который осуществляет конвекцию, B, A1, A2 — коэффициенты, которые необходимо идентифицировать. Для упрощения заменим все составляющие уравнения, которые мы можем считать константными в установившемся процессе и будем рассматривать следующее уравнение
![](https://habrastorage.org/r/w1560/getpro/habr/post_images/2a6/fdb/336/2a6fdb336ab4e90bb026072b69d59015.png)
Теперь у нас в распоряжении модель, от которой можно оттолкнуться при дальнейшем анализе данных. В эксперименте в качестве сигнала подавали меандр.
![](https://habrastorage.org/r/w1560/getpro/habr/post_images/afc/7a2/fd5/afc7a2fd57658eed53561818140e30b9.png)
![](https://habrastorage.org/r/w1560/getpro/habr/post_images/388/07a/366/38807a366361ab94e0da6d06937bf966.png)
Если мы получим оценку скорости температуры, то можно будет составить регрессионную модель с известными температурой и мощностью управляющего индуктора, и посчитать коэффициенты в уравнении.
Я уже достаточно давно использую для работы Scilab и в этот раз решил не изменять себе.Для этой задачи я выбрал взятие производной в образе Фурье. Но для начала работы с реальными данными, нужно интерполировать измерения для получения равномерной оси времени.
Стоит отметить, что переменная «int_temp_diff» будет содержать в себе данные скорости, но выглядят они весьма не лицеприятно.
![](https://habrastorage.org/r/w1560/getpro/habr/post_images/851/550/d91/851550d912e05ca5ae4e90c58d1f241d.png)
Так что перейдем к работе с данными в образе Фурье. Нам нужно будет нарастить дополнительный хвост данных, чтобы на конце не было разрыва, зеркально отобразив график данных температуры.
![](https://habrastorage.org/r/w1560/getpro/habr/post_images/71e/112/b32/71e112b32fcbe42fe75d62dfcff39634.png)
Потом задаем ось частот и делаем дискретное преобразование Фурье.
![](https://habrastorage.org/r/w1560/getpro/habr/post_images/bc3/d2f/b31/bc3d2fb31e0cfcf3388ac6366884c7c3.png)
![](https://habrastorage.org/r/w1560/getpro/habr/post_images/404/63f/fbb/40463ffbb7b797dd3be72ea246a3984a.png)
На частоте 1 герц виден явный пик, это связано с тем, что подложкодержатель вращается с этой частотой. Другие частоты пробились из-за того что графит нагрет неравномерно, и при вращении температура может увеличиться и уменьшиться несколько раз.
Чтобы взять производную в образе Фурье нужно просто домножить на
. После этого делаем обратное преобразование Фурье.
![](https://habrastorage.org/r/w1560/getpro/habr/post_images/5e6/c28/869/5e6c288695443140c7e2a7aa75c21eff.png)
Это конечно немного лучше, чем то что возвращает функция взятие производной через кубические сплайны, но идентификация с этим выдает некачественные результаты. Придется прибегнуть к фильтрации сигнала.
По совету моего друга Юры, который занимается настройкой корабельных систем управления, было решено использовать окно Блэкмана-Харриса для вырезания лишних частот. Так как в процессе роста графит вращается с частотой 1 раз в секунду, все колебания с частотой выше этой нам не интересны. Поэтому стоит вырезать все колебания выше 1 [Гц]:
![](https://habrastorage.org/r/w1560/getpro/habr/post_images/0a6/665/8df/0a66658dfeb9a5f3737b8b826f9498c2.png)
Стало намного лучше. С этим приступим к вычислению коэффициентов регрессии.
Чтобы получить хорошие данные, необходимо предварительно нормировать вектора данных. Для свободной константы в уравнении нужно задать вектор заполненный единицами.
Теперь можно приступать к нахождению значений констант. Составляем матрицу, содержащую значения трех переменных и константы. При это делим каждый из векторов данных на его длину. Для получения коэффициентов нужно умножить псевдообратную матрицу на нормированный вектор данных скорости изменения температуры.
И после этого необходимо вернуться к оригинальным значениям вектора коэффициентов.
В результате получились следующие значения a0 = [ — 1.073D-12, — 0.0029096, 0.0004617, 2.0969723 ], что достаточно близко к результату наших зарубежных коллег, хотя у них используется ламповый нагрев.
Вот теперь можно приступать к проверке результатов на полученной модели. Воспользуемся для этого пакетом виртуального моделирования Xcos Scilab. Так как у нас есть дифференциальное уравнение
, описывающее термический процесс в реакторе, и все коэффициенты к нему, то соберем следующую схему.
![](https://habrastorage.org/r/w1560/getpro/habr/post_images/aec/52c/15b/aec52c15b9cdd3e3059f83d6b117a11d.png)
Блоком «From Workspace» подаем экспериментальные данные мощности, предварительно для этого блока создав структуру V. Задаем частоту сэмплирования и время начало моделирования в блоке «Clock». Длину выходной структуры данных задаем равной входной. Конечное время моделирования должно равняться времени эксперимента. Составленную модель запускаем в тексте программы и сравниваем графики эксперимента и моделирования.
![](https://habrastorage.org/r/w1560/getpro/habr/post_images/7a0/d86/ad2/7a0d86ad2b93f9764501579dbaa20875.png)
![](https://habrastorage.org/r/w1560/getpro/habr/post_images/9eb/482/70a/9eb48270a6a7d59195db4d58d77825a6.png)
Вот еще тоже самое, но для гармонического возмущения с изменяющейся частотой.
![](https://habrastorage.org/r/w1560/getpro/habr/post_images/854/e48/c73/854e48c73e2f4a8ad9cef5900e4b149d.png)
![](https://habrastorage.org/r/w1560/getpro/habr/post_images/d06/978/219/d06978219e6bf5b4beee5e7efb435162.png)
Огромное спасибо научному руководителю Arastas за неоценимую помощь и участие в подготовке материала. Благодарю Заварина Евгения за предоставление доступа к установке и программную реализацию всех экспериментов.
Физики из СПб ФТИ им. Иоффе занимаются выращиванием нитрид галлиевых полупроводниковых структур, которые обладают неплохими показателями скорости носителей заряда при переходе и большим коэффициентом теплопроводности. Процесс роста такой структуры проходит при температуре 1000 С (1273 К) и атмосферном давлении. Все происходит в специальной камере, находящейся в герметичной зоне. При выращивании структуры весь объем реактора и герметичной зоны заполняется азотом. В процессе роста структуры подложкодержатель вращается с частотой один раз в секунду. Такие операции относятся к быстрым термическим процессам, скорость изменения температуры в которых варьируется от нескольких единиц до сотен градусов в секунду.
![](https://habrastorage.org/getpro/habr/post_images/1b0/9cc/cb4/1b09cccb471fa6f01aa7430f3d6f428b.jpg)
Моей задачей было управление температурой графитового подложкодержателя при помощи индуктивного нагрева.
Технические характеристики установки выглядят следующим образом. Для измерения температуры используется лазерный пирометр, снимающий данные в центре графита. Частота съема информации 10 раз в секунду, шаг измерения 1 градус. Значение мощности передаваемой графиту полагается прямо пропорциональным мощности на индукторе. У генератора управляющего индуктором имеется цифровой выход, по которому передаются значения напряжения, тока и мощности.
Для начала нужно было настроить регулятор температуры, чтобы в процессе роста не было сильных колебаний. С этой задачей справились достаточно быстро, но наше решение давало качественный результат только на высоких температурах. Для других процессов требовалось менять коэффициенты.
Так как эта работа вписалась в тематику моей кандидатской, то захотелось написать какой-нибудь хитрый алгоритм управления на основе анализа модели термических процессов. Для начала ознакомился с законом Стефана — Больцмана, мощность излучения абсолютно чёрного тела прямо пропорциональна площади поверхности и четвёртой степени температуры. С учетом конвекции можно записать уравнение для термических процессов в точке съема температуры
![](https://habrastorage.org/getpro/habr/post_images/d68/0c5/fa8/d680c5fa8af50b1fa28108256851355b.png)
где T — температура, P — мощность, Tc — температура окружающих объектов, которые нагревают графит собственным излучением, Ta — температура газа возле точки съема информации, который осуществляет конвекцию, B, A1, A2 — коэффициенты, которые необходимо идентифицировать. Для упрощения заменим все составляющие уравнения, которые мы можем считать константными в установившемся процессе и будем рассматривать следующее уравнение
![](https://habrastorage.org/getpro/habr/post_images/2a6/fdb/336/2a6fdb336ab4e90bb026072b69d59015.png)
Теперь у нас в распоряжении модель, от которой можно оттолкнуться при дальнейшем анализе данных. В эксперименте в качестве сигнала подавали меандр.
![](https://habrastorage.org/getpro/habr/post_images/afc/7a2/fd5/afc7a2fd57658eed53561818140e30b9.png)
![](https://habrastorage.org/getpro/habr/post_images/388/07a/366/38807a366361ab94e0da6d06937bf966.png)
Если мы получим оценку скорости температуры, то можно будет составить регрессионную модель с известными температурой и мощностью управляющего индуктора, и посчитать коэффициенты в уравнении.
Я уже достаточно давно использую для работы 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» будет содержать в себе данные скорости, но выглядят они весьма не лицеприятно.
![](https://habrastorage.org/getpro/habr/post_images/851/550/d91/851550d912e05ca5ae4e90c58d1f241d.png)
Так что перейдем к работе с данными в образе Фурье. Нам нужно будет нарастить дополнительный хвост данных, чтобы на конце не было разрыва, зеркально отобразив график данных температуры.
![](https://habrastorage.org/getpro/habr/post_images/71e/112/b32/71e112b32fcbe42fe75d62dfcff39634.png)
Потом задаем ось частот и делаем дискретное преобразование Фурье.
![](https://habrastorage.org/getpro/habr/post_images/bc3/d2f/b31/bc3d2fb31e0cfcf3388ac6366884c7c3.png)
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
![](https://habrastorage.org/getpro/habr/post_images/404/63f/fbb/40463ffbb7b797dd3be72ea246a3984a.png)
На частоте 1 герц виден явный пик, это связано с тем, что подложкодержатель вращается с этой частотой. Другие частоты пробились из-за того что графит нагрет неравномерно, и при вращении температура может увеличиться и уменьшиться несколько раз.
Чтобы взять производную в образе Фурье нужно просто домножить на
![](https://habrastorage.org/getpro/habr/post_images/460/2c1/b6f/4602c1b6fefd468d4b7a5dfa9d4f2cbf.png)
omega = (2*%pi*%i*frequency);
tmp = sp.*omega;
out = real(fft(tmp,1));
speed_temp = out(1:N/2);
![](https://habrastorage.org/getpro/habr/post_images/5e6/c28/869/5e6c288695443140c7e2a7aa75c21eff.png)
Это конечно немного лучше, чем то что возвращает функция взятие производной через кубические сплайны, но идентификация с этим выдает некачественные результаты. Придется прибегнуть к фильтрации сигнала.
По совету моего друга Юры, который занимается настройкой корабельных систем управления, было решено использовать окно Блэкмана-Харриса для вырезания лишних частот. Так как в процессе роста графит вращается с частотой 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);
![](https://habrastorage.org/getpro/habr/post_images/0a6/665/8df/0a66658dfeb9a5f3737b8b826f9498c2.png)
Стало намного лучше. С этим приступим к вычислению коэффициентов регрессии.
Чтобы получить хорошие данные, необходимо предварительно нормировать вектора данных. Для свободной константы в уравнении нужно задать вектор заполненный единицами.
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. Так как у нас есть дифференциальное уравнение
![](https://habrastorage.org/getpro/habr/post_images/a43/e64/9d3/a43e649d369cd2b8384725daf9c7658f.png)
![](https://habrastorage.org/getpro/habr/post_images/aec/52c/15b/aec52c15b9cdd3e3059f83d6b117a11d.png)
Блоком «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--");
![](https://habrastorage.org/getpro/habr/post_images/7a0/d86/ad2/7a0d86ad2b93f9764501579dbaa20875.png)
![](https://habrastorage.org/getpro/habr/post_images/9eb/482/70a/9eb48270a6a7d59195db4d58d77825a6.png)
Вот еще тоже самое, но для гармонического возмущения с изменяющейся частотой.
![](https://habrastorage.org/getpro/habr/post_images/854/e48/c73/854e48c73e2f4a8ad9cef5900e4b149d.png)
![](https://habrastorage.org/getpro/habr/post_images/d06/978/219/d06978219e6bf5b4beee5e7efb435162.png)
Огромное спасибо научному руководителю Arastas за неоценимую помощь и участие в подготовке материала. Благодарю Заварина Евгения за предоставление доступа к установке и программную реализацию всех экспериментов.