用于传感器的多元回归分析-微分分析

目的

测试微分分析的方法是否可以提高传感器的响应速度

简单分析

多元回归-微分分析/导数回归

关联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);

未完待续