Введение

Контур измерения температуры

Измерительная подсистема

#include <OneWire.h>
// OneWire DS18B20 Temperature Example
// http://www.pjrc.com/teensy/td_libs_OneWire.html
// The DallasTemperature library can do all this work for you!
// http://milesburton.com/Dallas_Temperature_Control_Library
//Вывод Data преобразователя DS18B20 подключаем к 7-му выводу //платформы Arduino
OneWire  ds(7);  
void setup(void) {
  Serial.begin(9600);
}
void loop(void) {
  byte i;
  byte present = 0;
  byte type_s = 0;
  byte data[12];
  byte addr[8];
  float celsius;
  ds.search(addr);  
  ds.reset();
  ds.select(addr);
  ds.write(0x44, 1);  
  delay(800);    
  present = ds.reset();
  ds.select(addr);
  ds.write(0xBE);   // Read Scratchpad
  for ( i = 0; i < 9; i++) {  
    data[i] = ds.read();
   }
  int16_t raw = (data[1] << 8) | data[0];
  celsius = (float)raw / 16.0;
  Serial.println(celsius);
}

Python интерфейс для обработки “кривой разгона” и её графической реализации

1. Управление и опрос виртуального COM порта через USB порт; здесь использованы функции библиотеки pySerial.
2. Интерактивное взаимодействие с оператором; здесь вводим необходимое количество измерений и номер последовательного порта (смотрим в Диспетчере Устройств ОС MS Windows). Если номер порта введен верно и не занят другими программами, на консоль и в графическое окно выводятся текущие измеряемые значения технологического параметра. По окончании измерений, на консоли выводится продолжительность времени измерений, период опроса датчика и указание – в каком файле txt находится таблица с результатами измерений, а в графическое окно выводится график измерений, дата, время, номер эксперимента, период опроса датчика.
3. Автоматическая регистрация результатов измерений; здесь на диске в папке со скриптом необходимо наличие файла count.txt, в котором записано целое без знаковое число, используемое как счётчик экспериментов; каждый последующий эксперимент увеличивает счётчик на единицу и добавляет значение счётчика к имени файла с таблицей результатов эксперимента.

import numpy as np
import matplotlib.pyplot as plt
import serial
from drawnow import drawnow
import datetime, time
#вывод выборки в графическое окно
def cur_graf():
    plt.title("DS18B20")
    plt.ylim( 20, 40 )
    plt.plot(nw, lw1, "r.-")
    plt.ylabel(r'$температура, \degree С$')
    plt.xlabel(r'$номер \ измерения$')
    plt.grid(True)
#вывод всех списков в графическое окно
def all_graf():
    plt.close()
    plt.figure()
    plt.title("DS18B20\n" +
              str(count_v) + "-й эксперимент " +
              "(" + now.strftime("%d-%m-%Y %H:%M") + ")")
    plt.plot( n, l1,  "r-")
    plt.ylabel(r'$температура, \degree С$')
    plt.xlabel(r'$номер \ измерения$' +
               '; (период опроса датчика: {:.6f}, c)'.format(Ts))
    plt.grid(True)
    plt.show()
#определяем количество измерений
# общее количество измерений
str_m   = input("введите количество измерений: ") 
m = eval(str_m)
# количество элементов выборки
mw  = 16
#настроить параметры последовательного порта
ser = serial.Serial()
ser.baudrate = 9600
port_num = input("введите номер последовательного порта: ")
ser.port = 'COM' + port_num
ser
#открыть последовательный порт
try:
    ser.open()
    ser.is_open
    print("соединились с: " + ser.portstr)
except serial.SerialException:
    print("нет соединения с портом: " + ser.portstr)
    raise SystemExit(1)
#определяем списки
l1  = [] # для значений температуры
t1  = []
lw1 = [] # для значений выборки температуры
n   = [] # для значений моментов времени
nw  = [] # для значений выборки моментов времени
#подготовить файлы на диске для записи
filename = 'count.txt'
in_file = open(filename,"r")
count = in_file.read()
count_v = eval(count) + 1
in_file.close()
in_file = open(filename,"w")
count = str(count_v)
in_file.write(count)
in_file.close()
filename = count + '_' + filename
out_file = open(filename,"w")
#вывод информации для оператора на консоль
print("\n параметры:\n")
print("n - номер измерения;")
print("T - температура, град. С;")
print("\n измеряемые значения величины температуры\n")
print('{0}{1}\n'.format('n'.rjust(4),'T'.rjust(10)))
#считывание данных из последовательного порта
#накопление списков
#формирование текущей выборки
#вывод значений текущей выборки в графическое окно
i = 0
while i < m:
    n.append(i) 
    nw.append(n[i])
    if i >= mw:
        nw.pop(0)
    line1 = ser.readline().decode('utf-8')[:-2]
    t1.append(time.time())
    if line1:
        l1.append(eval(line1))
        lw1.append(l1[i])
        if i >= mw:
            lw1.pop(0) 
    print('{0:4d} {1:10.2f}'.format(n[i],l1[i]))
    drawnow(cur_graf)
    i += 1
#закрыть последовательный порт
ser.close()
ser.is_open
#time_tm -= time_t0
time_tm = t1[m - 1] - t1[0]
print("\n продолжительность времени измерений: {0:.3f}, c".format(time_tm))
Ts = time_tm / (m - 1)
print("\n период опроса датчика: {0:.6f}, c".format(Ts))
#запись таблицы в файл
print("\n таблица находится в файле {}\n".format(filename))
for i in np.arange(0,len(n),1):
    count = str(n[i]) + "\t" + str(l1[i]) + "\n"
    out_file.write(count)
#закрыть файл с таблицей
out_file.close()
out_file.closed
#получить дату и время
now = datetime.datetime.now()
#вывести график в графическое окно
all_graf()
end = input("\n нажмите Ctrl-C, чтобы выйти ")

Python интерфейс для получения передаточной функции и оценки адекватности модели

# -*- coding: utf8 -*-
import matplotlib.pyplot as plt
import time
start = time.time()
from scipy.interpolate import splev, splrep
import scipy.integrate as spint
import numpy as np
from scipy.integrate import quad
xx =np.array(np.arange(0,230,1))
yy1 =np.array([
 24.87, 25.06, 25.31, 25.5, 25.81, 26.06, 26.37, 26.62, 26.87, 27.12, 27.44, 27.69, 27.87, 28.12, 28.31, 28.56, 28.75, 28.94, 29.12, 29.31, 29.5, 29.69, 29.81, 30.0, 30.12, 30.25, 30.37, 30.56, 30.69, 30.81, 30.87, 31.0, 31.12, 31.25, 31.31, 31.44, 31.5, 31.62, 31.69, 31.81, 31.87, 31.94, 32.06, 32.13, 32.19, 32.25, 32.31, 32.38, 32.44, 32.5, 32.56, 32.63, 32.69, 32.75, 32.75, 32.81, 32.88, 32.94, 32.94, 33.0, 33.06, 33.13, 33.13, 33.19, 33.25, 33.25, 33.31, 33.31, 33.38, 33.38, 33.44, 33.44, 33.5, 33.5, 33.56, 33.56, 33.56, 33.63, 33.63, 33.69, 33.69, 33.69, 33.75, 33.75, 33.81, 33.81, 33.81, 33.88, 33.88, 33.88, 33.88, 33.94, 33.94, 33.94, 34.0, 34.0, 34.0, 34.0, 34.06, 34.06, 34.06, 34.06, 34.06, 34.13, 34.13, 34.13, 34.13, 34.13, 34.19, 34.19, 34.19, 34.19, 34.19, 34.25, 34.25, 34.25, 34.25, 34.25, 34.25, 34.31, 34.31, 34.31, 34.31, 34.31, 34.31, 34.31, 34.38, 34.38, 34.38, 34.38, 34.38, 34.38, 34.38, 34.38, 34.38, 34.44, 34.44, 34.44, 34.44, 34.44, 34.44, 34.44, 34.44, 34.44, 34.5, 34.5, 34.5, 34.5, 34.5, 34.5, 34.5, 34.5, 34.5, 34.5, 34.5, 34.5, 34.5, 34.56, 34.56, 34.56, 34.56, 34.56, 34.56, 34.56, 34.56, 34.56, 34.56, 34.56, 34.56, 34.56, 34.56, 34.56, 34.56, 34.63, 34.56, 34.63, 34.63, 34.63, 34.63, 34.63, 34.63, 34.63, 34.63, 34.63, 34.63, 34.63, 34.63, 34.63, 34.63, 34.63, 34.63, 34.63, 34.63, 34.63, 34.69, 34.63, 34.69, 34.69, 34.69, 34.69, 34.69, 34.69, 34.69, 34.69, 34.69, 34.69, 34.69, 34.69, 34.69, 34.69, 34.69, 34.69, 34.69, 34.69, 34.69, 34.69, 34.69, 34.69, 34.69, 34.69, 34.75, 34.69, 34.75, 34.69, 34.69, 34.75, 34.75, 34.75, 34.75, 34.75])
yy2=yy1-24.87#компенсация смещения нуля
yy=yy2/max(yy2)#нормирование
""" Интерполяция переходной характеристики при помощи сплайнов"""
def h(x):         
         spl = splrep(xx , yy )
         return splev(x, spl)
""" Численное интегрирование без смены координаты времени в соответствии с (6)"""
S1=(spint.quad(lambda x:1-h(x),xx[0],xx[len(xx)-1])[0])
S2=(spint.quad(lambda x:(1-h(x))*(S1-x),xx[0],xx[len(xx)-1])[0])
S3=(spint.quad(lambda x:(1-h(x))*(S2-S1*x+(1/2)*x**2),xx[0],xx[len(xx)-1])[0])
S4=(spint.quad(lambda x:(1-h(x))*(S3-S2*x+S1*(1/2)*x**2-(1/6)*x**3),xx[0],xx[len(xx)-1])[0])
""" Определение коэффициентов передаточной функции"""
b1=-S4/S3
a1=b1+S1
a2=b1*S1+S2
a3=b1*S2+S3
""" Возврат во временную область"""
def ff(x,t):
        j=(-1)**0.5
        return (2/np.pi)*( ((b1*x*j+1)*np.e**(-25*x)/(a3*(x*j)**3+a2*(x*j)**2+a1*x*j+1)).real)*(np.sin(x*t)/x)
y=np.array([round(quad(lambda x: ff(x,t),0, 0.6)[0],2) for t in xx])
""" Определение критерия адекватности модели """
k=round(1-sum([(yy[i]-y[i])**2 for i in np.arange(0,len(yy)-1,1)])/sum([(yy[i])**2 for i in np.arange(0,len(yy)-1,1)]),5)
stop = time.time()
print ("Время работы программы :",round(stop-start,3))
plt.title('Идентификация модифицированным методом площадей.\n Адекватность модели: %s'%k)
plt.plot(xx, yy,label='W=(%s*p+1)/(%s*p**3+%s*p**2+%s*p+1)'%(round(b1,1),round(a3,1),round(a2,1),round(a1,1)))
plt.legend(loc='best')
plt.grid(True)
plt.show()

Выводы