Робототехника
0,0
рейтинг
6 апреля 2014 в 00:25

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

Недавно мне удалось завершить часть работы по очень интересному проекту в ФТИ им. Иоффе и получить достаточное количество экспериментальных данных, для того чтобы поделиться с Вами.
Физики из СПб ФТИ им. Иоффе занимаются выращиванием нитрид галлиевых полупроводниковых структур, которые обладают неплохими показателями скорости носителей заряда при переходе и большим коэффициентом теплопроводности. Процесс роста такой структуры проходит при температуре 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 за неоценимую помощь и участие в подготовке материала. Благодарю Заварина Евгения за предоставление доступа к установке и программную реализацию всех экспериментов.
Robot Engineering Department @kap2fox
карма
36,0
рейтинг 0,0
Робототехника
Реклама помогает поддерживать и развивать наши сервисы

Подробнее
Спецпроект

Самое читаемое Разработка

Комментарии (20)

  • +3
    Нет графиков разности (по модулю) между экспериментальными и полученными опытными данными с течением времени…
  • +4
    Спасибо за замечание, добавил.
  • +2
    ФТИ им. Йоффе

    Неужели так теперь пишут? Вроде был Иоффе
    • +1
      Уже исправил, спасибо)
  • +6
    Мне кажется, давно пора статью на хабре приравнять к ВАКовской публикации. Можно отдельный хаб для этого создать, а еще ввести специальную «научную» карму :)

    • +6
      Это востребованная и правильная идея, только нужно:
      1) Создать специальную группу экспертов по данной тематике. Только им позволить плюсовать/минусовать научную карму.
      2) Ввести порог научной кармы, после которой статья считается ВАКовской.
      • +1
        Это верно, но самым сложным будет добиться собственно решения ВАКа. Хотя, наверное, это не невозможно. Может быть у кого-нибудь есть там контакты, подбросить идею? Ну и инициатива, конечно, должна исходить от Тематических Медиа.
        • +1
          Я эту идею уже пару лет вынашиваю. Но в моей области (радиотехника) это пока нереально — не сидят наши корифеи в интернетах. А вот в области Computer Science экспертов на хабре найдется достаточно.

          Добиться решения ВАКа, на мой взгляд, не проблема:
          1) Насколько я помню, есть журналы, издаваемые электронно, но входящие в список ВАК. Группа экспертов хабра де-юре будет рецензентами. И, я уверен, куда более эффективными, чем рецензенты в печатных изданиях.
          2) Хабр в почете у наших представителей власти. А обсуждаемое — «инновационная» идея, так что поддержку будет легко найти.
        • +2
          Может быть сперва развить инфраструктуру, а потом уже начинать контакт с ВАК?
    • 0
      Вы знаете, мне кажется, что это два очень разных формата, в чем смысл их смешивать? Если есть хороший материал, то можно и пост сделать, и статью подготовить для тематического ВАКовского журнала.
      • 0
        Проблема в том, что ВАКовские журналы почти никто не читает. Получается работа ради «галочки». Нет feedback'а, нет обсуждения.
        • +1
          Предположим, ВАК станет учитывать некоторые топики на хабре как публикации. Что изменится? Начнут больше читать журналы ВАК? Люди (ученые), который сейчас не знают про хабр и не умеют пользоваться поисковиком, вдруг научатся и придут на хабр? В чем смысл?

          И потом, что значит — не читают? ВАКовские журналы бывают очень разные. Какой-то вестник/сборник трудов малоизвестного ВУЗа, наверное, не читают. Переводные журналы, входящие в WoS или Scopus — читают. Публикуйтесь во вторых или сразу в серьезных зарубежных, их больше читают. Хотите обсуждения — ездите на конференции, общайтесь, рассказывайте о своих результатах.
          • 0
            Это не должен быть хабр. Нужна некоторая электронная площадка, на которой можно было бы публиковаться как можно чаще. Без заморочек с рецензированием и т.д. Это позволило бы увеличить число статей, расширить круг вопросов, создало бы площадку для обсуждения, увеличило количество вовлеченных людей.

            Периодически бы проскакивали серьезные статьи, которые по мнению «старейшин» могли бы быть удостоены «ВАКовского» звания.

            Я работаю в области спутниковой радионавигации. Наши ВАКовские журналы выходят раз-два в год, публикации в них — скорее для отчетности. Иностранные журналы по нашей тематике тоже не радуют. Есть единственная нормальная конференция, американский ION, и их сборник статей. Но уходить играть на чужое поле — это не лучший вариант. Нужно чтобы и у нас была жизнь, а не те 5-10 статей из России в год, что сейчас публикуют американцы.

            Итого:
            1) Увеличить количество публикаций
            2) Увеличить скорость публикаций
            3) Получать обратную связь и обсуждение
            4) Вовлечь больше специалистов
            5) Иметь возможность интерпретировать топик как полноценную статью при соответствующем его уровне
            • +1
              Хотите много и часто публиковаться — пишите на arxiv.org. Бороться нужно за ужесточение рецензирования и предпечатной редактуры/корректуры, а не за их отмену. Собственно, наличие независимого рецензирования — один из главных признаков площадки для научных публикаций. По крайне мере в моей области (системы управления).
  • 0
    Нам нужно будет нарастить дополнительный хвост данных, чтобы на конце не было разрыва, зеркально отобразив график данных температуры.


    А зачем? Это ничего не дает. Для борьбы с разрывом надо не наращивать график, а домножить его на подходящую взвешивающую функцию («окно»). Отражение данных — это прием который используется для анализа функций, про которые из других источников известно что они четные или нечетные. Но Вы вообще просто дифференцированием занимаетесь, там совсем необходимости в дополнительных манипуляциях нет.
  • 0
    Емнип, применение оконных функций само по себе лишних частот не вырезает — оно изменяет влияние факта конечности длины блока данных на результат ПФ. Если же нужно вырезать все, что выше определенной частоты, можно рассчитать и применить цифровой ФНЧ нужного порядка
    • +1
      Автор, как я понимаю, окно применил к частотному спектру, а не исходному сигналу. Довольно странная идея, кстати — специально несовершенный фильтр НЧ получился. А ФНЧ напрямую не гонялся потому как автору его хотелось к производной применить, а для этого производную нужно было вначале посчитать, что делалось работой с Фурье-образом.

      Впрочем для выбранного фильтра что производную вначале посчитать, а затем ФНЧ прогнать, что в обратной последовательности ФНЧ прогнать, а затем производную считать — разницы не будет.
  • 0
    Не думали вместо преобразования Фурье с оконной функцией использовать вейвлет-преобразование, например с базисом Добеши?
    • 0
      Сейчас Мы уже получили необходимые оценки математической модели, которые дают хороший результат. В будущем обязательно попробуем предлженный метод. Спасибо за совет.
  • 0
    На мой взгляд, для оценивания скорости удачным выбором является предварительная фильтрация, а затем аппроксимация экспериментальной кривой сглаживающими сплайнами и аналитическое нахождение производной у аппроксимирующей кривой. Недостатком будет необходимость руками выбирать параметры сглаживания при аппроксимации и, возможно, вручную разделить исходные данные на несколько интервалов с разным сглаживанием.
    Фильтрацию же в оффлайне можно сделать и через Фурье, и через двойное прохождение ФНЧ в прямом-обратном времени.

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