目的
测试微分分析的方法是否可以提高传感器的响应速度
简单分析
多元回归-微分分析/导数回归
关联6个变量
浓度值 C | EMF1 | EMF1' | EMF2 | EMF2' | 时间 t |
输出浓度值=常数项+AE1+BE1'+CE2+DE2'+平均误差
A~D为试验数据计算得到的回归系数
单片机ADC采样得到的E1、E2为离散的变量,采用前向差分来求得导数的近似值
Python 代码实现
简单测试
import itertools
# 计算转置矩阵
def transpose(matrix):
return list(map(list, zip(*matrix)))
# 计算矩阵乘法
def matrix_multiply(A, B):
return [[sum(a * b for a, b in zip(A_row, B_col)) for B_col in zip(*B)] for A_row in A]
# 计算矩阵的逆
def matrix_inverse(matrix):
size = len(matrix)
I = [[float(i == j) for i in range(size)] for j in range(size)]
for i in range(size):
factor = matrix[i][i]
if factor == 0:
for k in range(i + 1, size):
if matrix[k][i] != 0:
matrix[i], matrix[k] = matrix[k], matrix[i]
I[i], I[k] = I[k], I[i]
factor = matrix[i][i]
break
else:
raise ValueError("矩阵不可逆")
for j in range(size):
matrix[i][j] /= factor
I[i][j] /= factor
for k in range(size):
if k != i:
factor = matrix[k][i]
for j in range(size):
matrix[k][j] -= factor * matrix[i][j]
I[k][j] -= factor * I[i][j]
return I
# 测试数据
t = [1,6,12,18,24,30,36,42,48,54]
E1 = [944,901,851,804,759,726,689,665,640,613]
E2 = [664,670,665,664,662,658,655,655,653,648]
C = [10,20,30,40,50,60,70,80,90,100]
# 计算导数
def forward_difference(data, time):
derivatives = []
for i in range(len(data) - 1):
derivative = (data[i + 1] - data[i]) / (time[i + 1] - time[i])
derivatives.append(derivative)
last_derivative = (data[-1] - data[-2]) / (time[-1] - time[-2])
derivatives.append(last_derivative)
return derivatives
E1_prime = forward_difference(E1, t)
E2_prime = forward_difference(E2, t)
# 将数据整合到一个矩阵中
X = []
for i in range(len(t)):
X.append([1, E1[i], E1_prime[i], E2[i], E2_prime[i]]) # 包含常数项1
Y = C
# 转换为矩阵形式
X_matrix = X
Y_matrix = [[y] for y in Y]
# 计算 (X^T * X)^(-1) * X^T * Y
X_transpose = transpose(X_matrix)
X_transpose_X = matrix_multiply(X_transpose, X_matrix)
X_transpose_X_inv = matrix_inverse(X_transpose_X)
X_transpose_Y = matrix_multiply(X_transpose, Y_matrix)
beta = matrix_multiply(X_transpose_X_inv, X_transpose_Y)
# 展开结果
beta = [b[0] for b in beta]
# 打印回归系数
print("regression coefficient:", beta)
print("Constant term:", beta[0])
feature_names = ["E1", "E1'", "E2", "E2'"]
for idx, coef in enumerate(beta[1:], start=1):
print(f"{feature_names[idx-1]} coefficient:", coef)
# 计算误差项(残差)
Y_pred = [sum(b * x for b, x in zip(beta, row)) for row in X_matrix]
residuals = [y - y_pred for y, y_pred in zip(Y, Y_pred)]
print("Error term:", residuals)
输出
regression coefficient: [1441.4044757131633, -0.1288821345484097, 0.541399530490537, -1.9528115291560155, -7.377922099411137]
Constant term: 1441.4044757131633
E1 coefficient: -0.1288821345484097
E1' coefficient: 0.541399530490537
E2 coefficient: -1.9528115291560155
E2' coefficient: -7.377922099411137
Error term: [0.43665714164164804, 1.465446387396117, -0.09480301810906511, 0.48480493556112236, -2.762620604880759, -3.2363904580402334, -1.3475352901566566, 3.190219369292791, 2.5540484077285868, -0.6898268708586244]
Avr Error term: -4.2507330988428293e-11
测试ok
后续测试流程
应用到单片机
假设差分求导时间间隔为1s
关键问题在于求E1E2关于时间的导数
单片机实现
if (voltage_EMF1<6)//待修改
{
if (timer_3min>=10&&Previous_voltage==false)
{
V_EMF_1=voltage_EMF1*10000*0.1;
V_EMF_2=voltage_EMF2*10000*0.1;
PRE_V_EMF_1=V_EMF_1;
PRE_V_EMF_2=V_EMF_2;
Previous_voltage=true;
}
if (timer_3min>=20&&Previous_voltage==true)
{
D_V_EMF_1=(PRE_V_EMF_1-V_EMF_1)/0.1;
D_V_EMF_2=(PRE_V_EMF_2-V_EMF_2)/0.1;
timer_3min=0;
Previous_voltage=false;
}
else if (timer_3min>30)
{
timer_3min=0;
Previous_voltage=false;
}
concentration=(1398.4441-0.5241*V_EMF_1-0.1187*D_V_EMF_1-1.4966*V_EMF_2-0.0766*D_V_EMF_2+0)*10;
if (concentration<0)
{
concentration=0;
}
}
else
{
concentration=0;
}
Matlab实现测试数据转换为计算式系数
实际应用需要很多数据导入
求导间隔=1s
% 读取数据
data = readtable('C:\Users\Leon\Desktop\6.xlsx'); % 文件路径
% 提取各列数据
t = data{:, 1}; % 第1列为时间t
E1 = data{:, 2}; % 第2列为E1
E2 = data{:, 3}; % 第3列为E2
C = data{:, 4}; % 第4列为响应变量C
% 计算导数(差分法)
E1_prime = diff(E1) ./ diff(t);
E2_prime = diff(E2) ./ diff(t);
% 为导数补上最后一个元素,使其长度与原数据一致
E1_prime = [E1_prime; E1_prime(end)];
E2_prime = [E2_prime; E2_prime(end)];
% 构造回归矩阵X,包含常数项1
X = [ones(length(t), 1), E1, E1_prime, E2, E2_prime];
% 转换响应变量为列向量
Y = C;
% 计算回归系数 (X^T * X)^(-1) * X^T * Y
beta = (X' * X) \ (X' * Y);
% 打印回归系数
fprintf('回归系数:\n');
fprintf('常数项: %f\n', beta(1));
feature_names = {'E1', 'E1''', 'E2', 'E2'''};
for idx = 2:length(beta)
fprintf('%s 系数: %f\n', feature_names{idx-1}, beta(idx));
end
% 计算预测值
Y_pred = X * beta;
% 计算残差
residuals = Y - Y_pred;
% 打印残差
fprintf('残差:\n');
disp(residuals');
% 计算平均残差
mean_residual = mean(residuals);
fprintf('平均残差: %f\n', mean_residual);
求导间隔=1s
% 读取数据
data = readtable('C:\Users\Leon\Desktop\6.xlsx'); % 文件路径
% 提取各列数据
t = data{:, 1}; % 第1列为时间t,时间间隔为0.1s
E1 = data{:, 2}; % 第2列为E1
E2 = data{:, 3}; % 第3列为E2
C = data{:, 4}; % 第4列为响应变量C
% 计算导数(差分法),时间间隔为0.1s
time_interval = 0.1;
E1_prime = diff(E1) / time_interval;
E2_prime = diff(E2) / time_interval;
% 为导数补上最后一个元素,使其长度与原数据一致
E1_prime = [E1_prime; E1_prime(end)];
E2_prime = [E2_prime; E2_prime(end)];
% 构造回归矩阵X,包含常数项1
X = [ones(length(t), 1), E1, E1_prime, E2, E2_prime];
% 转换响应变量为列向量
Y = C;
% 计算回归系数 (X^T * X)^(-1) * X^T * Y
beta = (X' * X) \ (X' * Y);
% 打印回归系数
fprintf('回归系数:\n');
fprintf('常数项: %f\n', beta(1));
feature_names = {'E1', 'E1''', 'E2', 'E2'''};
for idx = 2:length(beta)
fprintf('%s 系数: %f\n', feature_names{idx-1}, beta(idx));
end
% 计算预测值
Y_pred = X * beta;
% 计算残差
residuals = Y - Y_pred;
% 打印残差
fprintf('残差:\n');
disp(residuals');
% 计算平均残差
mean_residual = mean(residuals);
fprintf('平均残差: %f\n', mean_residual);
未完待续