Поиск

Кусочно-линейная аппроксимация


Вчера я занимался построением зависимости отчетов в ADU от четырех каналов измеряемой температуры контроллера чиллера. Для этого замотал на термопасте четыре NTC и один платиновый терморезистор (для контроля) между двух кусочков фольгированного стеклотекстолита. Поместил это в баночку из-под детского питания, залил баночку тосолом, поставил в банку из-под кофе и залил азотом. Как только азот испарился (температура тосола опустилась до -140°C), залил в железную банку теплого тосола и стал фиксировать значения ADU и температуры по мере нагревания.

В итоге у меня получилась табличка, которую теперь нужно обработать. Каждая запись в таблице — медианное усреднение девяти последовательных измерений с интервалом в 0.1с. В среднем СКО получилось около 2.5 на все четыре канала. На низких температурах расхождение их показаний не сильно отличалось от СКО, однако, по мере уменьшения сопротивления расхождение было все больше и больше. Несмотря на то, что используемые в качестве опорных резисторов имеют паспортную точность в 0.1% (т.е. 1Ом), их сопротивления имеют значение от 994 до 998 Ом, вот, собственно, где собака и порылась!
Я не нашел, как в octave найти узлы кусочно-линейной интерполяции, поэтому пришлось велосипедить.
Итак, я занес данные в файл Temperatures.dat, где первый столбец — сопротивление терморезистора, второй — вычисленная по этому сопротивлению температура, а оставшиеся четыре — показания АЦП контроллера чиллера. Так как мне не нужна температура ниже -20°C, первые строки можно обрезать:

 d = dlmread('Temperatures.dat'); d(1:17,:) = []; T = d(:,2); % температуры R = d(:,[3:6]); % ADU сопротивлений


Строим матрицу для кубической аппроксимации

 M3=[ones(size(T)) T T.^2 T.^3];


И вычисляем коэффициенты:

 km3 = M3\median(R,2) km3 = 1002.2286743 37.1381302 0.3308417 -0.0044131


Строим новые матрицы с шагом 0.05 градусов:

 Tt = [-20:0.05:20]'; Mt = [ones(size(Tt)) Tt Tt.^2 Tt.^3]; ADU = Mt * km3;


Теперь имеем обратную зависимость Tt(ADU), которой нужно подобрать кусочно-линейную аппроксимацию.

Ищем узлы:

 knots = [1 getknots(ADU, Tt, 0.1) max(size(ADU))] knots = 1 42 84 169 257 350 446 550 664 801

Строим график:

 plot(ADU, Tt, ADU(knots), Tt(knots), '+') print -dpng knots.png

Формируем таблицу:

 [X0 Y0 K]=buildk(ADU, Tt, 0.1) X0 = 427.11 467.72 514.28 622.83 753.63 909.75 1087.41 1295.45 1537.77 Y0 = -20.0000 -17.9500 -15.8500 -11.6000 -7.2000 -2.5500 2.2500 7.4500 13.1500 K = 0.050477 0.045107 0.039150 0.033639 0.029785 0.027017 0.024996 0.023522 0.022514

Сверяем:

 mR2 = median(R,2); cT=[];for i = 1:max(size(mR2)); cT=[cT calcT(mR2(i))]; endfor plot(R,T,'o', mR2,T,'*', mR2,cT); print -dpng check.png

В области высоких температур все получается довольно-таки хреново:

размах показаний такой, что вытягивает почти на 1°C! Однако, найденная кусочно-линейная интерполяция вполне дает желаемую точность в ±0.5°C. Таким образом, строить индивидуальные зависимости для каждого канала не нужно! Кстати, если ограничиться размахом в 0.5°C (а не в 0.1°C), то получится всего четыре узла! Но, т.к. 9 узлов не так-то много места в памяти займут, пусть будет так!
Смеха ради посмотрел, сколько будет узлов с размахом не больше 0.01°C. Всего-то 28 штук! Т.е. данная методика вполне подходит и для аппроксимации измерений температур при помощи платиновых терморезисторов на внешнем АЦП. Можно будет ею воспользоваться, если доберусь-таки до системы точной регистрации температур (узел АЦП/мультиплексора я еще летом спаял, нужно лишь подключить терморезисторы и микроконтроллер, да написать прошивку).

Исходники утилит.

function Tout = H705(Rin)
% Converts resistance of TRD into T (degrC) _alpha = 0.00375;
_beta = 0.16;
_delta = 1.605;
T = [-300:0.1:300];
_A = _alpha + _alpha*_delta/100.;
_B = -_alpha*_delta/1e4;
_C = zeros(size(T));
_C(find(T<0.)) = -_alpha*_beta/1e8;
rT = 1000.*(1 + _A*T + _B*T.^2 - _C.*T.^3*100. + _C.*T.^4);
%plot(T, rT);
Tout = interp1(rT, T, Rin, 'spline');
endfunction

function idx = getknots(X, Y, dYmax)
%
% idx = getknots(X, Y, dYmax) - calculate piecewise-linear approximation knots for Y(X)
% dYmax - maximal deviation % idx = []; I = getnewpt(X, Y, dYmax); if(I > 0) L = getknots(X(1:I), Y(1:I), dYmax); R = I-1 + getknots(X(I:end), Y(I:end), dYmax); idx = [L I R]; endif
endfunction

function idx = getnewpt(X, Y, delt)
%
% idx = getnewpt(X, Y, delt)
% find point where Y-`linear approx` is maximal
% return index of max deviation or -1 if it is less than `delt`
% % make approximation [x0 y0 k] = linearapprox(X,Y); % find new knot adiff = abs(Y - (y0 + k*(X-x0))); maxadiff = max(adiff); if(maxadiff < delt) idx = -1; else idx = find(adiff == max(adiff)); endif
endfunction

function [X0 Y0 K] = buildk(X, Y, dYmax)
%
% [X0 Y0 K] = buildk(X, Y, dYmax) - build arrays of knots (X0, Y0) and linear koefficients K
% for piecewise-linear approximation of Y(X) with max deviation `dYmax`
% X0 = []; Y0 = []; K = []; knots = [1 getknots(X, Y, dYmax) max(size(X))]; for i = 1:max(size(knots))-1 idx = knots(i):knots(i+1); [x y k] = linearapprox(X(idx), Y(idx)); X0 = [X0 x]; Y0 = [Y0 y]; K = [K k]; endfor
endfunction

function [x0 y0 k] = linearapprox(X,Y)
% % [a b] = linearapprox(X,Y) - find linear approximation of function Y(X) through first and last points
% y = y0 + k*(X-x0)
% x0 = X(1); y0 = Y(1); k = (Y(end) - y0) / (X(end) - x0);
endfunction

И функция для проверки (можно было бы построить структуру при помощи mkpp и дальше средствами octave работать, но мне интересно было проверить наиближайшее приближение к тому, что будет на МК):

function T = calcT(ADU)
%
% T = calcT(ADU) - piecewise calculation prototype
% X0 = [427.11 467.72 514.28 622.83 753.63 909.75 1087.41 1295.45 1537.77]; Y0 = [-20.0000 -17.9500 -15.8500 -11.6000 -7.2000 -2.5500 2.2500 7.4500 13.1500]; K = [0.050477 0.045107 0.039150 0.033639 0.029785 0.027017 0.024996 0.023522 0.022514]; imax = max(size(X0)); idx = int32((imax+1)/2); T = []; % find index while(idx > 0 && idx <= imax) x = X0(idx); half = int32(idx/2); if(ADU < x) %printf("%g < %g\n", ADU, x); if(idx == 1) break; endif; % less than first index if(ADU > X0(idx-1)) idx -= 1; break; endif; % found idx = half; % another half else %printf("%g > %g\n", ADU, x); if(idx == imax) break; endif; % more than last index if(ADU < X0(idx+1)) break; endif; % found idx += half; endif endwhile if(idx < 1) printf("too low (<%g)!\n", X0(1)); idx = 1; endif if(idx > imax) idx = imax; endif; T = Y0(idx) + K(idx) * (ADU - X0(idx));
endfunction


eddy_em.livejournal.com

Добавить комментарий