METHOD FOR THE EARLY ESTIMATION OF ANAEROBIC DEGRADABILITY OF ORGANIC SUBSTRATES

Abstract

A method for the early estimation of anaerobic degradability of organic substrates, starting from initial data acquired from tests measuring BMP (Biochemical Methane Potential). The method consists of: i) calculating the two parameters B.sub.0 and k; ii) comparing the fit of the decreasing trend of B.sub.0,est as Δt varies with a homographic function in the first quadrant; iii) evaluating the goodness of fit between a homographic function and the trend of B.sub.0,est as Δt varies, checking whether the adjusted coefficient of determination R.sup.2.sub.adj≥R.sup.2.sub.adj,min; iv) selecting the value of B.sub.0,est corresponding to a slope of less than 0.1% that occurs for three consecutive Δt; if no, acquire additional BMP measurements and repeat the previous steps.

Claims

1. A method for the early estimation of anaerobic degradability of organic substrates comprising the following steps: performing a periodic sampling of the BMP measurement; acquiring a pre-set number of BMP measurements; considering that the anaerobic degradability of organic substrates follows the following formula:
Y(t)=B.sub.0(1−e.sup.−k(t−t.sup.0.sup.)) where Y(t) is the value of the signal at time t, B.sub.0 is the maximum production of methane at infinite time, k is the kinetic constant of the first order; t.sub.0 is the time at which acquisition of the measurements begins; calculating the two parameters B.sub.0 and k, measured at time t.sub.1, beginning of measurements, and at time (t.sub.1+Δt), where Δt is a generic time interval and is identified as B.sub.0,est and k.sub.est; the time interval Δt is incremented, and calculation of B.sub.0,est and k.sub.est is repeated until Δt reaches a pre-set value Δ.sub.tmax; comparing the fit of the decreasing trend of the curve B.sub.0,est as Δt varies with a homographic function in the first quadrant; evaluating whether the fit between a homographic function and trend of the curve B.sub.0,est as Δt varies, has an adjusted coefficient of determination R.sup.2ad.sub.j greater than or equal to a pre-set minimum coefficient of determination R.sup.2.sub.adj,min; analysing the curve of B.sub.0,est as Δt varies, and select the value of B.sub.0,est corresponding to a slope of the curve of B.sub.0,est of less than 0.1% that occurs for three consecutive Δt; if no, acquire additional BMP measurements and repeat the previous steps.

2. The method in accordance with claim 1 characterised in that it performs a moving average with a window W at said acquired data.

3. The method in accordance with claim 1 characterised in that the step of calculating the two parameters B.sub.0 and k takes place by means of the method described in the article by Niemann & Miklos, which calculates the parameters B.sub.0 and k, once a part of the initial response of the system y(t) and its prime derivative y′(t) are known, given by the following equations:
y(t)=B.sub.0.Math.(1−e.sup.−k.Math.(t−t.sup.0.sup.))
y′(t)=k.Math.B.sub.0.Math.e.sup.−k.Math.(t−t.sup.0.sup.) where t.sub.0 is the time instant in which the measurements are started.

4. The method in accordance with claim 3 characterised in that the step of calculating the two parameters B.sub.0 and k comprises the step of calculating a first prime derivative of a spline calculated in first three predefined points.

5. The method in accordance with claim 3 characterised in that the step of calculating the two parameters B.sub.0 and k comprises the step of calculating a second prime derivative of a spline calculated in second three predefined points.

6. The method in accordance with claim 4 characterised in that the step of calculating the two parameters B.sub.0 and k comprises the step of calculating the average between said first derivative and said second derivative.

7. The method in accordance with claim 1 characterised in that said homographic function is a generic function having the following equation: $y=\frac{a.Math.x+b}{c.Math.x+d}$ where the parameters a, b, c, d are estimated using nonlinear least squares regression techniques.

8. The method in accordance with claim 1 characterised in that said adjusted coefficient of determination R.sup.2.sub.adj is given by the following equation: $R_{\mathrm{adj}}^2=\left[\frac{\left(1-R^2\right).Math.\left(n-1\right)}{n-z-1}\right]$ where R.sup.2 is the coefficient of determination, n is the number of data available, and z is the number of independent variables of the regression.

9. The method in accordance with claim 1 characterised in that the step of selecting the value of B.sub.0,est corresponding to a slope of less than 0.1% that occurs for three consecutive zit comprises evaluating the following conditions:
[B.sub.0,est(i)−B.sub.0,est(i−1)]/B.sub.0,est(i)≤slope.sub.max
[B.sub.0,est(i−1)−B.sub.0,est(i−2)]/B.sub.0,est(i−1)≤slope.sub.max
[B.sub.0,est(i−2)−B.sub.0,est(i−3)]/B.sub.0,est(i−2)≤slope.sub.max where index i identifies the progressive number of estimations of B.sub.0,est that were calculated.

10. The method in accordance with claim 1 characterised in that the step of calculating the two parameters B.sub.0 and k takes place for given values of t.sub.1 and Δt, and where the generic time interval Δt is initially initialized at t.sub.1.

Description

[0016] The features and advantages of the present invention will become clear from the following detailed description of a practical embodiment thereof, illustrated by way of non-limiting example in the accompanying drawings, in which:

[0017] FIG. 1 shows a flowchart of a method for the early estimation of anaerobic degradability of organic substrates, in accordance with the present invention;

[0018] FIG. 2, in support of FIG. 1, shows how the two parameters t.sub.min and p, which are used for the early estimation of anaerobic degradability of organic substrates, are defined through an example provided for a generic data vector having a length equal to 30, in accordance with the present invention;

[0019] FIG. 3 shows graphs representing the calculation as set forth in the article by Niemann & Miklos, in accordance with the present invention.

[0020] In accordance with the present invention, the method is based on the fact that, in a large number of real-world cases, the kinetics of substrate conversion to methane predominantly follows dynamics of the first order as can be seen in the following formula:

(t)=B.sub.0.Math.(1−e.sup.−k.Math.(t−t.sup.0.sup.))

[0021] where y(t) is the value of the signal at time t, t.sub.0 is the test start time, B.sub.0 is the maximum methane production at infinite time and k is the kinetic constant.

[0022] First order kinetics includes the cumulative effects of all the biochemical reactions that take place during the anaerobic degradation process, assuming the initial hydrolysis process as limiting. For this reason, and given its simplicity, which implies a limited number of parameters to be estimated, the dynamics of the first order remains the most widely used model to describe BMP curves.

[0023] In addition, a mathematical method proposed in the following article will be used to estimate the two parameters B.sub.0 and k, from which it is possible to calculate the value of BMP at the end of the test.

[0024] Niemann, H. H., & Miklos, R. (2014). A Simple Method for Estimation of Parameters in First order Systems. Journal of Physics: Conference Series (Online), 570, [012001]. https://doi.org/10.1088/1742-6596/570/1/012001.

[0025] Assuming that a generic number m of test data is acquired, once a generic estimation time t.sub.1 and a generic time interval Δt (such that Δt≤m−t.sub.1) are defined, the procedure described by Niemann & Miklos estimates the first order kinetic parameters by using only part of the experimental data, and presuming that the prime derivative of the function at the time instant t.sub.1 is known. This method, if applied to BMP tests, would allow to obtain very low estimation errors for specific values of t.sub.1 and Δt, which, however, are not known a priori and must be appropriately selected to minimize the errors in the parameter estimations. In fact, due to the presence of noise in the experimentally acquired data, and dynamics that are not always attributable to a pure model of the first order, it is possible that, as the pairs t.sub.1 and Δt vary, the more or less significant estimation errors of the parameters B.sub.0 and k alternate.

[0026] In the method developed it has been attempted to optimize the aspects considered most critical such as the removal of noise from the raw data acquired experimentally; the reliable estimation of the prime derivative of the curve at t.sub.1; the identification of a method to select a given combination t.sub.1−Δt at which the parameters B.sub.0 and k are estimated early with an error that is considered acceptable.

[0027] The method developed allows, given a certain number of data acquired on an hourly basis from an experimental test, to select the pair of values t.sub.1−Δt at which to estimate a final BMP value (derived from the parameters B.sub.0 and k) which has a good probability of reducing the estimation error within values considered as acceptable. Through the method indicated, the kinetic parameters of the first order model are progressively calculated for each possible pair t.sub.1−Δt: when specific conditions are reached in the trend of the estimations of B.sub.0 as Δt varies, the method stops the iterations and returns the estimated BMP value.

[0028] The method comprises the following steps.

[0029] 1). Acquire 12 and store a number of acquisitions every hour, n.sub.data, of the measurements performed by an ongoing BMP test 10. A minimum number of data is required to start the algorithm (n.sub.data,min). This number n.sub.data,min is defined, equal to at least 25, including the initial zero instant.

[0030] 2). When the minimum number n.sub.data,min of measurements has been acquired, the condition that n.sub.data≥n.sub.data,min is checked 14. If yes Y, the calculation method starts, and the number of data h.sub.tot subsequently used for calculations is initialized 16 to n.sub.data,min (h.sub.tot=n.sub.data,min).

[0031] If no N, the acquisitions 12 of the performed BMP measurements 10 continue, and the number n.sub.data is incremented by one unit, step by step, until the pre-set total number of data n.sub.data,min is reached.

[0032] When n.sub.data,min is reached, storing stops (but data acquisition continues) and the variables B.sub.0 and k (denoted by B.sub.0,est and k.sub.est where the suffix est stands for estimated value) are estimated early as will be seen in the next steps; h.sub.tot also represents the actual number of data used to obtain the estimation, correlated with the estimation time necessary for parameter prediction.

[0033] 3). In order to remove the noise present in the acquired data 18 a moving average with a window W of length preferably equal to 7 samples (equivalent to 7 hours), is applied to the data. The new data vector will therefore have the same length as the original (h.sub.tot); however, at the beginning and at the end of this vector, the reprocessed data cannot be considered reliable as the window W of the moving average is unbalanced at the extremes. For this reason, it is necessary to define two time indices in order to select the corresponding BMP data values, y(t), excluding the extremes of this vector.

[0034] A minimum value train, equal to 7, defined in order to properly calculate and average the value y(t.sub.int,1), where the index t.sub.int,1 is used to calculate the prime derivative at t.sub.1, as will be further explained below.

[0035] A maximum index, t.sub.max, calculated as the difference between h.sub.tot and p, where p is equal to 3 and takes into consideration the number of data that cannot be appropriately averaged when approaching the upper extreme of the vector of available data. It follows that with p=3, the last 7 values are averaged correctly.

[0036] 4). It is therefore now possible to start calculating, for each possible pair t.sub.1−Δt, the early estimations of B.sub.0 and k. In order to reduce the number of iterations to what is actually necessary with respect to subsequent processing, the maximum time index that can be assumed by t.sub.1 (equal to t.sub.1,max) is defined 20 as the difference between t.sub.max and v(B.sub.0).sub.min, where v(B.sub.00.sub.min=25; it follows that t.sub.1,max is a function of only t.sub.max since v(B.sub.0).sub.min always remains constant, the latter being the minimum number of early estimations of B.sub.0 necessary to decide whether or not to stop the iterations and then select the estimated parameters.

[0037] 5). The generic index t.sub.1 is initialized 22 at t.sub.min.

[0038] 6). Given the value of t.sub.1, the maximum time interval Δt.sub.max to be considered is defined 24 as the difference between t.sub.max and t.sub.1.

[0039] 7). Subsequently, the generic time interval Δt is initialized 26 at t.sub.1, and B.sub.0,est and k.sub.est are calculated 28 based on the experimental values of methane produced measured at t.sub.1 and at (t.sub.1+Δt), and based on the value of the prime derivative of the experimental curve calculated at t.sub.1. The calculation is performed as indicated in the article by Niemann & Miklos.

[0040] The method, allows the estimation of the first order kinetic parameters (k and B.sub.0), once a part of the initial response of the system y(t) and its prime derivative y′(t) are known, given by the following equations:

y(t)=B.sub.0.Math.(1−e.sup.−k(t−t.sup.0.sup.))

y′(t)=k.Math.B.sub.0.Math.e.sup.−k(t−t.sup.0.sup.)

[0041] where t.sub.0 is the time instant in which the step function or initial step is applied to the system. The curves of y(t) and y′(t) can be seen in the graph on the left-hand side of FIG. 3.

[0042] Once a generic time instant t.sub.1 and a generic time interval Δt are defined, and assuming they are known: [0043] the values of the response signal of the system at t.sub.1 and at t.sub.1+Δt, defined as y(t.sub.1) and y(t.sub.1+Δt), respectively; [0044] the value of the prime derivative at t.sub.1, defined as y′(t.sub.1);

[0045] it is possible to derive the following equation that allows to evaluate the constant k by solving the aforementioned equation for k:

[00001]$R_1\left(\Delta t\right)=\frac{y\left(t_1+\Delta t\right)-y\left(t_1\right)}{y^{\prime }\left(t_1\right).Math.\Delta t}=\frac{1-e^{-k.Math.\Delta t}}{k.Math.\Delta t}=\frac{1-e^{-\alpha }}{\alpha }=R\left(\alpha \right)$

[0046] Therefore, the calculation of R.sub.1(Δt) is performed with the data found in the graph on the left-hand side of FIG. 3, and k is represented by the graph on the right-hand side of FIG. 3, where α is in turn the product between the kinetic constant k (unknown) and the time interval Δt adopted (known). As the value of R.sub.1(Δt) is known, it is in fact univocally possible to estimate the value of α, by solving the following equation for α:

R.sub.1(Δt)−R(α)=0

[0047] Instead, the constant Bo is calculated as shown in the following equation:

B.sub.0=y(t.sub.1)+y′(t.sub.1)/k

[0048] The prime derivative of the experimental curve calculated at t.sub.1, is the result of the average between two values of prime derivative at t.sub.1, calculated: in the first case starting from a prime derivative at t.sub.1 of a spline evaluated in the points t.sub.int,1− t.sub.1−t.sub.int,2; in the second case, starting from a prime derivative at ti of a spline evaluated in the points t.sub.1−t.sub.int,2−t. Where t is equal to the sum of t.sub.1 and Δt, and t.sub.int,1 and t.sub.int,2 are the intermediate points between t.sub.0−t.sub.1 and t.sub.1−t, respectively.

[0049] 8). The time interval Δt is then incremented by one time unit (one hour) 30, and calculation 28 of B.sub.0,est and k.sub.est is repeated 32 until Δt reaches Δt.sub.max.

[0050] 9). To stop the calculation and select the parameters B.sub.0 and k, it is necessary to check 34 the fitting to a homographic function in the first quadrant, of the decreasing trend of the curve of B.sub.0,est as Δt varies.

[0051] A homographic function is a generic function having the following equation:

[00002]$y=\frac{a.Math.x+b}{c.Math.x+d}$

[0052] where the parameters a, b, c, d are estimated by nonlinear least squares regression techniques. In the present case, the Levenberg-Marquardt algorithm was used.

[0053] It has in fact been observed that, at a generic value of t.sub.1, the trend of the curve of B.sub.0,est as a function of Δt is similar to that of the estimation error of the parameter itself, and has a decreasing trend as explained. On the other hand, the estimations of k and their errors do not always have the same trend or at least an easily identifiable trend: for this reason, selecting the correct estimation of k can be more complex.

[0054] 10). It is checked whether the fit 36 is good between homographic function and trend of B.sub.0,est as Δt varies, i.e., R.sup.2adj R.sup.2adj,min with R.sup.2.sub.adj,min=0,998 and where R.sup.2.sub.adj is the adjusted coefficient of determination, namely:

[00003]$R_{\mathrm{adj}}^2=\left[\frac{\left(1-R^2\right).Math.\left(n-1\right)}{n-z-1}\right]$

[0055] where R.sup.2 is the coefficient of determination, n is the number of data available, and z is the number of independent variables of the regression.

[0056] If the fit 36 is not good, i.e., R.sup.2.sub.adj<R.sup.2.sub.adj,min, t.sub.1 is incremented 42 by one unit and until 44 t.sub.1≤t.sub.1,max the operations from 24 to 34 are performed.

[0057] It should be specified that the generic time interval Δt is initialized at t.sub.1 in order to remove the initial noise due to very small values of Δt.

[0058] 11). It is checked 38 whether the function of B.sub.0,est has a slope of less than 0.1%, defined as slope.sub.max=0.001, i.e., that slope≤slope.sub.max of the function B.sub.0,est that occurs for three consecutive Δt.

[0059] That is, the following conditions must hold simultaneously:

[B.sub.0,est(i)−B.sub.0,est(i−1)]/B.sub.0,est(i)≤slope.sub.max

[B.sub.0,est(i−1)−B.sub.0,est(i−2)]/B.sub.0,est(i−1)≤slope.sub.max

[B.sub.0,est(i−2)−B.sub.0,est(i−3)]/B.sub.0,est(i−2)≤slope.sub.max

[0060] where index i identifies the progressive number of estimations of B.sub.0,est that have been calculated.

[0061] If yes, B.sub.0,est is selected 40 (as is the corresponding value of k.sub.est) and the method stops and returns the result.

[0062] Otherwise, the value of tris incremented by one 42 and the operations from step 24 to step 34 are repeated.

[0063] 12). If it is still not possible to select any value of B.sub.0,est and k.sub.est, and 44 (N) t.sub.1>t.sub.1,max, h.sub.tot is incremented 46 by one and, if 48 h.sub.tot≤n.sub.data, step 18 and the subsequent steps are then performed, which are repeated until all the data acquired by the instrument are used (n.sub.data), i.e., until 48 h.sub.tot is less than or equal to n.sub.data. If no (N), and if no solution is found, this means that more data 10 are required before a reliable early estimation of the constants B.sub.0 and k and therefore the BMP can be obtained.

[0064] In accordance with the present invention, it should be specified that the time necessary for the estimation does not depend upon the kinetic constant k; the time necessary for the estimation, and the goodness and reliability of the same, do not depend upon the type of substrate tested, nor upon the magnitude of the final BMP value found; the estimation errors of B.sub.0 are less than or at least comparable with what has already been shown in the literature.

[0065] The method, as it is structured, can be easily integrated into a data acquisition software for BMP tests.

[0066] The method has been developed, and contextually optimized and checked, on multiple BMP test cases. To make this study more robust, note that these tests were conducted using different testing setup parameters (I/S ratio) and inoculum biomasses having different origins. In addition, the substrates tested are grouped into four different categories: primary sewage sludge, secondary sewage sludge, agricultural waste/scrap, foodstuff waste/scrap. Through the use of Akaike and Bayesian criteria, it was finally verified that the prevailing dynamics of the experimental curve of each of the BMP cases selected was a first order kinetics.

[0067] The method thus conceived is susceptible to numerous modifications and variations, all falling within the inventive concept; furthermore, all the details can be replaced by technically equivalent elements.