Effect of Humidity on Electrochemical Sensor Signal and Seeking Correction Approaches
Show older comments
Hi everyone,
I've uploaded my dataset, which consists of a timetable with a model predictions, REF values (which my model should ideally match closely), humidity, and temperature readings. I've also included the plots I've generated based on this data.
The dataset I'm working with represents corrected signals from an electrochemical sensor. I've applied regression modeling and temperature correction to this data. Additionally, for the purpose of this analysis, I've removed outliers to ensure the data's reliability. Upon examining the plots, it's evident that as humidity levels fluctuate, they influence the overall signal. I believe there's a correlation here, especially since the impact spans over several hours, necessitating hours of data monitoring for accurate prediction.
While the sensor's output closely follows the trend of the reference signal, noticeable discrepancies arise in line with humidity changes. Specifically, the signal tends to overestimate when humidity increases and underestimate as it decreases. I've attempted to address this using the Regression Learner App and by comparing the error with humidity readings. However, I've not yet been successful in deriving a suitable correction factor or model. The challenge seems to be the need to assess prolonged data periods and discern if the humidity is rising or falling during those times.
Could anyone provide guidance or suggest an approach to tackle this issue?
Accepted Answer
Relative humidity will of course be affected by temperature, inbcreases as the temperature decreases, and it can never exceed 100%. It might be better to express humidity in terms of the dewpoint (the environmental temperature can never fall below the dewpoint) or absolute humidity, although this is also a function of temperature.
Plotting the data (‘Temperature’, ‘Humidity’, and the difference between ‘Predicted_DATA’ and ‘REF’) fails (to me,a t least) to reveal any patterns I might use to model them.
LD = load('REG_Model_Data.mat');
T1 = LD.cleaned_timetable_REG
VN = T1.Properties.VariableNames;
delta = T1{:,3}-T1{:,2};
Lv = delta >= 0;
figure
scatter3(T1{Lv,4}, T1{Lv,5}, delta(Lv), '.g')
hold on
scatter3(T1{~Lv,4}, T1{~Lv,5}, delta(~Lv), '.r')
patch([xlim flip(xlim)], [[1 1]*min(ylim) [1 1]*max(ylim)], zeros(1,4), [1 1 1]*0.5, 'FaceAlpha',0.25)
hold off
grid on
xlabel(VN{4})
ylabel(VN{5})
zlabel([VN{3} ' - ' strrep(VN{2},'_','\_')])
view(-60,30)
figure
scatter(T1.Temperature, T1.Humidity,'.')
grid
xlabel('Temperature (°C)')
ylabel('Relative Humidity (%)')
% delta = TT1.Predicted_DATA - TT1.REF;
% pred = [TT1.Humidity TT1.Temperature];
% mdl = fitlm(pred, delta)
I would include humidity as an additional independent variable in your model (if you have not already) and proceed from there. I doubt that there is any way to predict it, other than noting the rate and direction of change in both humidity and temperature.
.
23 Comments
Dharmesh Joshi
on 24 Oct 2023
Star Strider
on 24 Oct 2023
Edited: Star Strider
on 25 Oct 2023
That most likely has to do with the physical dynamics of the sensor. It might be possible to model the humidity changes in real time, and then have the sensor do minor corrections based on that, however it is not obvious to me how to incorporate that otherwise. Perhaps using a state-space representation of the system, where inputs are the instantaneous humidity and NOx concentrations, and the output is a correct NOx concentration could work. The state-space model would take into consideration the changing dynamics and correct for them.
Beyond that, I have no suggestions other than that the System Identification Toolbox could be helpful in estimating that system from the known inputs and known output.
EDIT — (25 Oct 2023 at 01:22)
I gave some thought to the state space model of temperature and humidity to predict the ‘REF’ data and came up with this preliminary result. The model is large for such models (this one is order 8, and a smaller system might do reasonably well — your decision, not mine) however if you are interested in this sort of approach, this should get you started.
The idea here is to account for the dynamics of the sytem as a function of time, not simply the instantaneous values of the input variables. Since you want to track the changes in humidity, similar approaches (such as a Kalman filter) might also be appropriate to estimate the directions of the various inputs to your system. (I do not have sufficient experience with Kalman filters or knowledge of your system and instrumentation to attempt that.)
The inputs here are ‘Humidity’ and ‘Temperature’ (in that order), and the output is.‘REF’.
A preliminary approach —
LD = load('REG_Model_Data.mat');
T1 = LD.cleaned_timetable_REG
figure
plot3(T1.Humidity, T1.Temperature, T1.REF)
grid on
xlabel('Humidity')
ylabel('Temperature')
zlabel('REF')
figure
plot(T1.Time, T1{:,[4 5 3]})
grid
legend('Humidity','Temperature','REF', 'Location','best')
SSData = iddata(T1{:,3}, T1{:,[4 5]}, 1, 'TimeUnit','hours')
sys = ssest(SSData, 8)
figure
compare(SSData, sys)
grid
.
Dharmesh Joshi
on 25 Oct 2023
Thank you; I will review the suggested documents.
In terms of humidity, its effects differ from temperature as they are not immediate but may take a few hours to manifest. This delay depends on whether the humidity is increasing or decreasing and at what intensity. Consequently, incorporating humidity into the standard model doesn't significantly impact, based on my experience. Visually comparing the rise and fall of humidity over a span of several hours reveals a correlation. To illustrate this more clearly, I have attached a screenshot of the data collected over a few days.
I am not certain that this analysis adds much, however it might be of interest —
LD = load('REG_Model_Data.mat');
T1 = LD.cleaned_timetable_REG
mdl = fitlm(T1.Humidity, T1.REF)
figure
plot(mdl) % x1 = Humidity, y = REF
grid
[r, lags] = xcorr(T1.Humidity, T1.REF, 'coeff'); % Return Correlation Coefficient
[rmax,idx] = max(r);
Delay = lags(idx)
figure
stairs(lags, r)
grid
xlabel('Hours')
ylabel('R_{x,y}')
text(lags(idx), rmax, sprintf('\\leftarrow R_{x,y} = %.3f\n \\Delta = %d hours', rmax, Delay), 'Horiz','left', 'Vert','top')
There does appear to be a correlation, and no significant delay.
If these are not the variables you want to examine, substitute the correct variables.
.
Dharmesh Joshi
on 25 Oct 2023
Thanks i will give that a go.
For your information this what is mentioned on the datasheet. But as humdity is change constantly, this transient issue is much longer.
Star Strider
on 25 Oct 2023
My pleasure.
The state space approach should have a way to deal with the transients, however it needs to have them as part of its data in order to model them effectively. That (or some version of it) may do what you want.
Dharmesh Joshi
on 30 Oct 2023
Hi
I have cleaned up my original signal (per minute data sample) as mentioned in my previous post.MATLAB Tools to Evaluate and Mitigate Sensor Noise - MATLAB Answers - MATLAB Central (mathworks.com)
For your reference, I have attached both the figure and raw data. The raw data comprises the REF data (hourly), raw data (per minute), and temperature/humidity (hourly). The goal is to correct the raw data so that, in principle, it looks similar to the REF signal; I am not expecting it to match 100%. I am currently evaluating the original 1-minute signal before it undergoes any type of regression modeling. Previously, I had provided the hourly averaged signal, but it was altered by a regression model. This time, it's the raw signal that I aim to correct due to humidity. Please see the plots below. In the plot, the blue represents REF data, black represents the raw data, and cyan represents humidity.
As you can observe, the raw signal decreases as humidity increases, and vice versa: it increases as humidity decreases. The intensity of these changes also correlates with the fluctuations in the raw signal.
To address this issue, I am exploring methods to correct the signal using a compensation function rather than creating a regression model at this stage. What approach would you recommend for developing such a compensation function? Ideally using the REF as abousulte signal, but use it as reference shaping signal
I belive after this correction, we could so a simple linear regression to the REF data for final value.
Star Strider
on 30 Oct 2023
I do not see any new files. Are they supposed to be here, or simply the plot images? (I just need to know.)
Dharmesh Joshi
on 30 Oct 2023
I have been looking at this for a while (interrupted by computer problems), and note that ‘smooth_difference’ lags ‘REF’ by 15 hours, and leads ‘RH’ by 9.
In estimating this using the System Identification Toolbox (which is the only approach I consider appropriate here), what should the input (or inputs) be, and what should the output be? I do not completely understand what you want to do in that respect.
I used retime to resample (interpolate) ‘smooth_diffference’ to the time vectors shared by all the other signals (I checked to be sure that they all share the same time, so I could use it with the other signals), and created a new structure for it
My analysis so far —
LD = load('plotData_needs...orrection.mat')
plot1 = LD.plot1
plot2 = LD.plot2
plot3 = LD.plot3
plot4 = LD.plot4
P2 = timetable(plot2.Time, plot2.smooth_difference); % Resample 'plot2' To Common 'Time'
P2r = retime(P2, plot1.Time);
plot2r.Time = P2r.Time;
plot2r.smooth_difference = P2r{:,1} % Create New Structure
figure
hold on
plot(plot1.Time, plot1.RH, 'DisplayName',plot1.DisplayName)
plot(plot2r.Time, plot2r.smooth_difference, 'DisplayName',plot2.DisplayName)
plot(plot3.Time, plot3.REF, 'DisplayName',plot3.DisplayName)
plot(plot4.Time, plot4.T, 'DisplayName',plot4.DisplayName)
hold off
grid
legend('Location','northoutside', 'NumColumns',2)
D1 = finddelay(plot1.RH, plot2r.smooth_difference) % Units: Hours
D2 = finddelay(plot3.REF, plot2r.smooth_difference) % Units: Hours
D3 = finddelay(plot3.REF, plot1.RH) % Units: Hours
% maxlag = 24;
% Dn = finddelay(plot3.REF, [plot1.RH; plot2r.smooth_difference; plot4.T], maxlag)
figure
hold on
% plot(plot1.Time, plot1.RH, 'DisplayName',plot1.DisplayName)
plot(plot2r.Time, plot2r.smooth_difference*4, 'DisplayName',plot2.DisplayName)
plot(plot3.Time, plot3.REF, 'DisplayName',plot3.DisplayName)
% plot(plot4.Time, plot4.T, 'DisplayName',plot4.DisplayName)
hold off
grid
xlim(plot1.Time(1)+[0 caldays(14)])
legend('Location','northoutside', 'NumColumns',2)
I believe it to be appropriate to model this as a dynamic system in state space using system identification, based on my understanding of what you are doing (that may not be correct).
.
Dharmesh Joshi
on 31 Oct 2023
I am having a difficult time visualising this. It seems that the input to the system are all four variables collected simultaneously, and thge output is a correction signal that is (‘REF’ - ‘smooth_difference’).
Testing that hypothesis —
LD = load('plotData_needs...orrection.mat')
plot1 = LD.plot1
plot2 = LD.plot2
plot3 = LD.plot3
plot4 = LD.plot4
P2 = timetable(plot2.Time, plot2.smooth_difference); % Resample 'plot2' To Common 'Time'
P2r = retime(P2, plot1.Time);
plot2r.Time = P2r.Time;
plot2r.smooth_difference = P2r{:,1} % Create New Structure
figure
hold on
plot(plot1.Time, plot1.RH, 'DisplayName',plot1.DisplayName)
plot(plot2r.Time, plot2r.smooth_difference, 'DisplayName',plot2.DisplayName)
plot(plot3.Time, plot3.REF, 'DisplayName',plot3.DisplayName)
plot(plot4.Time, plot4.T, 'DisplayName',plot4.DisplayName)
hold off
grid
legend('Location','northoutside', 'NumColumns',2)
SysData = iddata((plot2r.smooth_difference-plot3.REF), [plot1.RH plot2r.smooth_difference plot3.REF plot4.T], 1, 'TimeUnit','hours')
SSsys = ssest(SysData, 12)
figure
compare(SysData, SSsys)
grid
Tsim = 1:height(plot1.Time);
correction_vector_output = lsim(SSsys, [plot1.RH plot2r.smooth_difference plot3.REF plot4.T], Tsim);
smooth_difference_corrected = plot3.REF + correction_vector_output;
figure
plot(plot2r.Time, plot2r.smooth_difference, 'DisplayName','Original')
hold on
plot(plot2r.Time, smooth_difference_corrected, 'DisplayName','Corrected')
plot(plot3.Time, plot3.REF, 'DisplayName','REF')
hold off
grid
legend('Location','best')
figure
plot(plot2r.Time, plot2r.smooth_difference, 'DisplayName','Original')
hold on
plot(plot2r.Time, smooth_difference_corrected, 'DisplayName','Corrected')
plot(plot3.Time, plot3.REF, 'DisplayName','REF')
hold off
grid
legend('Location','best')
xlim([0 caldays(14)]+plot3.Time(1))
This actually appears to calculate the correction signal ‘(smooth_difference - REF)’ reasonably well. To get the corrected value of ‘smooth_difference’ add the ‘REF’ vector to it —
syms sd REF
correction_vector_output = sd - REF
smooth_difference_corrected = REF + correction_vector_output
I was hoping that the corrected signal would come out closer to ‘REF’ than to the original signal. I may have missed something in your description, so experiment with this approach, since you understand the problem much better than I do.
.
Dharmesh Joshi
on 31 Oct 2023
Please see the plot below. This is the trend I am referring to.
This is a zoomed-in view.
You can observe that the reference signal (in blue) is fairly flat. Then, there's the humidity signal (in cyan), which rises and falls throughout the day, affecting the smooth_difference signal. You'll notice that the smooth_difference signal increases its amplitude during the periods when the humidity drops, and does the opposite when it rises.
This is the impact of humidity. I have also included temperature data in case it provides any additional insights, but temperature can be ignored for now.
Additionally, I will review the steps you outlined to see if they help address the issue.
Dharmesh Joshi
on 1 Nov 2023
Star Strider
on 1 Nov 2023
It is difficult for me to understandconsider how to approach this. I still believe that a dynamic system approach is appropriate, however I am just not certain how to construct it. I do not have a clear understanding of what the inputs and outputs should be, since I do not understand the system you want to model.
Dharmesh Joshi
on 1 Nov 2023
Star Strider
on 1 Nov 2023
I had to run this offline because it takes longer than the allotted 55 seconds (it needs about 4½ times that).
load('REG_Model_Data.mat')
cleaned_timetable_REG
% NrMissing = nnz(ismissing(cleaned_timetable_REG))
% cleaned_timetable_REG = fillmissing(cleaned_timetable_REG, 'nearest')
cleaned_timetable_REG_TT = table2timetable(cleaned_timetable_REG);
cleaned_timetable_REG_TT = retime(cleaned_timetable_REG_TT, 'hourly');
cleaned_timetable_REG_TTs = cleaned_timetable_REG_TT(:,[1 3 4 2])
Data = iddata(cleaned_timetable_REG_TTs)
Data = misdata(Data)
SSsys = ssest(Data, 8)
figure
compare(Data, SSsys)
grid
It produces this system:
SSsys =
Continuous-time identified state-space model:
dx/dt = A x(t) + B u(t) + K e(t)
y(t) = C x(t) + D u(t) + e(t)
A =
x1 x2 x3 x4 x5 x6 x7 x8
x1 -0.09591 0.1753 0.1319 0.1004 -0.07117 -0.1091 -0.01476 -0.1967
x2 0.09113 -0.1586 -0.2084 -0.09554 0.1261 0.1156 0.02182 0.2023
x3 -0.1483 0.2734 -0.06624 -0.4684 0.2209 0.154 0.02348 0.1978
x4 0.06936 -0.09895 0.4544 -0.08761 0.2892 0.1671 0.06083 0.2402
x5 0.04081 -0.06385 -0.0145 -0.02059 -0.178 -0.9404 0.07352 -1.255
x6 0.04111 -0.04807 -0.141 -0.1098 0.7128 -0.1866 -1.424 -0.09405
x7 0.01679 -0.01909 0.01431 -0.03315 -0.04821 1.374 0.005255 -1.766
x8 0.02773 -0.04175 -0.09524 -0.08299 0.7382 -0.4315 1.609 -0.3404
B =
Predicted_DA Humidity Temperature
x1 -4.047e-05 3.267e-05 0.0001242
x2 3.706e-05 -2.963e-05 -0.0001252
x3 4.425e-05 -8.677e-06 5.426e-05
x4 6.102e-05 -2.236e-05 5.346e-05
x5 2.801e-05 -1.11e-05 -2.314e-07
x6 -6.455e-06 3.922e-06 1.297e-05
x7 4.199e-06 -1.922e-06 -1.146e-06
x8 -8.82e-07 1.541e-06 1.093e-05
C =
x1 x2 x3 x4 x5 x6 x7 x8
REF -1751 160.9 -33.37 43.05 20.1 10.11 9.504 2.154
D =
Predicted_DA Humidity Temperature
REF 0 0 0
K =
REF
x1 -0.0005337
x2 0.0002469
x3 -0.0001914
x4 0.0002226
x5 0.0001757
x6 5.587e-06
x7 0.0001112
x8 4.167e-05
Parameterization:
FREE form (all coefficients in A, B, C free).
Feedthrough: none
Disturbance component: estimate
Number of free coefficients: 104
Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.
Status:
Estimated using SSEST on time domain data "Data".
Fit to estimation data: 87.19% (prediction focus)
FPE: 19.05, MSE: 18.78
Model Properties
Elapsed time is 231.098728 seconds.
The inputs are:
Predicted_DATA
Humidity
Temperature
and the output is:
REF
The order 8 system does a reasonably accurate approximation.
What would you want from a model such as this?
.
Dharmesh Joshi
on 1 Nov 2023
Star Strider
on 1 Nov 2023
It will output something approximating ‘REF’ since that is how I designed it. My impression is that is what you want, given those inputs. It does not output a compensation signal, although I attempted exactly that earlier, although with ‘REF’ also as an input.
I can experiment with the same sort of difference output however not with ‘REF’ as an input if that makes sense. The ‘Predicted_DATA’ would still be an input however, since that is what you specified. The output would be the difference signal vector, as previously. The only difference from the previous approach would be that ‘REF’ would not be an input.
Dharmesh Joshi
on 2 Nov 2023
Star Strider
on 2 Nov 2023
It can output a vector of differences, similar to the approach in my earlier Comment. That approach models the time-related dynamics of the system, including the transients.
If the ranges of all the variables are accounted for in the data available here, and if ‘REF’ will not always be available, then the regression approach might be the best option. The next step would be to determine the relationships abetween the variables. If you want to model them all as linear, then that would be straightforward (multiple linear regression) and would return coefficients for each variable, for example —
If the relationships were nonlinear (and that depends on the physics and chemistry involved, and could be a polynomial relation or a more complicated functional relation), that would have to be modeled for each of them and combined into a single function. If you have a mathematical description of these relationships, estimating the parameters would be straightforward. Those models and parameters could then be incorporated into a full model and given the inputs, would produce a single output that would (within a specific error tolerance that could also be estimated and modeled), be the result. This approach might not model the transients well, however it could model a steady-state system.
Dharmesh Joshi
on 2 Nov 2023
Star Strider
on 2 Nov 2023
It might also be possible to use the gradient function with each of the associated variables as an additional input (covariate). It computes the numeric derivative (in this instance with respect to time), however like all high-pass filters (and the gradient or derivative is one such), it accentuates the noise, so noise has to be minimised first by smoothing the signal. (I prefer using the sgolayfilt function for this, using a 3-order polynomial and varying the ‘framelen’ until I get the result I want. The smoothdata function is also an option.)
If the data are not regularly sampled or the sampling interval is not known, compute the derivative as:
dydt = gradient(y) ./ gradient(t);
(assuming vector arguments of the same dimensions) and proceed from there.
Plotting the variables and their derivatives as functions of time will demonstrate if this is a suitable option.
It might also be appropriate to calculate and plot the derivatives of the viariables with respect to each other. This might provide insight into their relationships, even if those results are not included in the model.
.
More Answers (0)
Categories
Find more on Correlation Models in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
