A Takagi–Sugeno-Kang Fuzzy Model Formalization of Eelgrass Leaf Biomass Allometry with Application to the Estimation of Average Biomass of Leaves in Shoots: Comparing the Reproducibility Strength of the Present Fuzzy and Related Crisp Proxies

Abstract

The identification of the functional relationship that regulates the variation of individual leaf biomass in terms of related area in eelgrass, allows the derivation of convenient proxies for a nondestructive estimation of the average biomass of the leaves in shoots. The concourse of these assessment methods is fundamental for assessing the performance of restoration efforts for this species that are based on transplanting techniques. Prior developments proposed proxies for a nondestructive estimation of aforementioned average biomass of leaves in shoots derived from allometric models for the dependence of leaf biomass in terms of linked area. The reproducibility power of these methods is highly dependent on analysis method and data quality. Indeed, previous results show that allometric proxies for average biomass of leaves in shoots produced by parameter estimates fitted from quality controlled data via nonlinear regression yield the highest reproducibility strength. Nevertheless, the use of data processing entails subtleties mainly related to the subjectivity of the criteria for the rejection of inconsistent replicates in raw data. Here we introduce efficient- data quality control- free surrogates derived from a first order Takagi-Sugeno-Kang fuzzy model aimed to approximate the mean response of eelgrass leaf biomass depending on associated area. A comparison of the performances of the allometric and the fuzzy model constructs identified using available raw data shows that the Takagi-Sugeno-Kang paradigm for individual leaf biomass in terms of related area produced the most precise proxies for observed average biomass of leaves in shoots. The present results show how gains derived from the outstanding approximation capabilities of the first order Takagi-Sugeno-Kang fuzzy model for the nonlinear dynamics can be extended to the realm of eelgrass allometry.

Appendices

Appendix 1: Data Quality Control Approach

Previous results demonstrate that it is reasonable to assume that for eelgrass the leaf area-to-leaf weight variation pattern should distribute about a power function-like trend [7]. Thus, by analyzing in a leaf area-to-weight plot, the spreading of leaf biomass values we can detect values that present severe deviations to a dominant power function trend. Since in the present settings leaf area is obtained by means of the length times width proxy [22] unduly deviations of replicates from an expected trend could be attributed to errors in leaf length or width estimations, imprecise gear for dry weight assessment or even due to inappropriate registering of measurements. In order to conceive a data cleaning criteria, we firstly observed that the set of $n_{rw}$ leaves constituting the raw data can be arranged into several groups $G\left( a \right) = \left\{ {w_{i} \left( a \right) | 1 \le i \le n\left( g \right)} \right\}$ formed by the set of $n\left( g \right)$ leaf biomass replicates $w_{i} \left( a \right)$ that associate to an observed leaf area value $a$. Considering that the median of a group of data is totally immune to the sample size and a robust estimator of scale, for groups of ten or more replicates, we adapted a Median Absolute Deviation (MAD) data cleaning procedure [7, 14]. For each one of the groups $G\left( a \right)$ acquired median is denoted by means of the symbol ${\text{MED}}\left\{ {w_{1} \left( a \right), \ldots ,w_{n\left( g \right)} \left( a \right)} \right\}$ or simply by means of ${\text{MED}}\left( {G\left( a \right)} \right)$ for short. Then, for each replicate $w_{j} \left( a \right)$ in $G\left( a \right)$ its absolute deviation $\delta_{j} \left( a \right)$ relative to the group median ${\text{MED}}\left( {G\left( a \right)} \right)$ is

$$\delta_{i} \left( a \right) = \left| {w_{i} \left( a \right) - {\text{MED}}\left( {G\left( a \right)} \right)} \right|$$

(27)

Similarly, we obtained the median of the set of absolute deviations denoted by the symbol $MED\left\{ {\delta_{1} \left( a \right), \ldots ,\delta_{n\left( g \right)} \left( a \right)} \right\}$. Following Huber [15] and by also eliciting that eelgrass leaf biomass values are log-normally distributed [7] then obtained Median Absolute Deviation of a group $G\left( a \right),$ denoted through ${\text{MAD}}\left( {G\left( a \right)} \right)$ is given by

$${\text{MAD}}\left( {G\left( a \right)} \right) = b{\text{MED}}\left\{ { \delta_{1} \left( a \right), \ldots ,\delta_{n\left( g \right)} \left( a \right)} \right\},$$

(28)

where $b = 1/Q\left( {0.75} \right)$, and $Q\left( {0.75} \right)$ stands for the $0.75$ quantile associated to a lognormal distribution. Henceforth, for the deletion of uneven replicates in a group $G\left( a \right)$ we used the decision criterion

$$\begin{aligned} {\text{MED}}\left( {G\left( a \right)} \right) & - T \cdot {\text{MAD}}\left( {G\left( a \right)} \right) < w_{j} \left( a \right) < \\ {\text{MED}}\left( {G\left( a \right)} \right) & + T \cdot {\text{MAD}}\left( {G\left( a \right)} \right), \\ \end{aligned}$$

(29)

where $w_{j} \left( a \right)$ stands for the $jth$ leaf in the group $G\left( a \right)$ and $T$ is the rejection threshold, which after Miller [23] we set at a value $T = 3.$

Appendix 2: Lin’s Concordance Correlation Coefficient $\left( \rho \right)$

The value of $\rho$ the concordance correlation coefficient [12, 24] is commonly used to assess how well a new set of observations $y$ mimic an original set $x$. In other words, the value of $\rho$ provides a measure of reproducibility as it evaluates the agreement between the variables $x$ and $y$ by assessing the extent to which they fall on the 45° line through the origin. Its value is characterized in terms of the ratio of the expected orthogonal squared distance from the diagonal $y = x$ to the expected orthogonal squared distance from the diagonal $y = x$ assuming independency. When $\rho$ is computed on a m-length data set (i.e., two vectors ($x_{1} ,x_{2} , \ldots ,x_{m}$) and ($y_{1} ,y_{2} , \ldots ,y_{m}$) the resulting statistics is denoted by means of $\hat{\rho }$ and calculated through

$$\hat{\rho } = \frac{{2s_{xy} }}{{s_{x}^{2} + s_{y}^{2} + \left( {\bar{x} + \bar{y}} \right)^{2} }},$$

(30)

being

$$\bar{x} = \frac{1}{m}\sum\limits_{j = 1}^{m} {x_{j} }$$

(31)

$$s_{x}^{2} = \frac{1}{m}\sum\limits_{j = 1}^{m} { (x_{j} - \bar{x})^{2} }$$

(32)

$$s_{y}^{2} = \frac{1}{m}\sum\limits_{j = 1}^{m} { (y_{j} - \bar{y})^{2} }$$

(33)

and

$$s_{xy} = \frac{1}{m}\sum\limits_{j = 1}^{m} {(x_{j} - \bar{x})(y_{j} - \bar{y})}$$

(34)

In the present work the value of $\hat{\rho }$ provides a criterion to assess to what extent the allometric proxies $w_{m} \left( {\alpha ,\beta ,t} \right)$ or the weighted output $w_{TSK} \left( {a(t)} \right)$ of the considered. Takagi–Sugeno-Kang model reproduce observed $w_{m} \left( t \right)$ values. Agreement will be defined as poor whenever $\hat{\rho } <$ 0.90, moderate for $0.90 \le \hat{\rho } <$ 0.95, good for $0.95 \le \hat{\rho } \le 0.99$, or excellent for $\hat{\rho } > 0.99$ [25].

Appendix 3: The Takagi-Sugeno-Kang Fuzzy Inference System

A Takagi-Sugeno-Kang inference system is a special characterization of what is known as a fuzzy inference system. Generally in setting a fuzzy inference system for $1 \le j \le n$ we initially consider input values $x_{j}$ an input domain $U$. The collection of input values is denoted by the symbol

$$X = \left\{ {x_{1} , \ldots ,x_{n} } \right\}$$

(35)

and similarly, for $1 \le s \le p$ we consider output values $y_{s} \left( {x_{1} , \ldots ,x_{n} } \right)$ in a range $V$, and the collection of output variables will be denoted by the symbol

$$Y = \left\{ {y_{1} \left( {x_{1} , \ldots ,x_{n} } \right), \ldots , y_{p} \left( {x_{1} , \ldots ,x_{n} } \right)} \right\}$$

(36)

For each input variable $x_{j}$ we associate a set $L_{j}$ containing a number $a\left( j \right)$ of linguistic terms $A_{k} \left( {x_{j} } \right)$, namely

$$L_{j} = \left\{ {A_{k} \left( {x_{j} } \right) |k = 1,2, \ldots ,q\left( j \right)} \right\}$$

(37)

The symbol $A_{X}$ will stand for the collection of linguistic terms that associate to $X$, formally

$$A_{X} = \cup_{j} L_{j} ,$$

(38)

The linguistic term $A_{k} \left( {x_{j} } \right)$ associates to a membership function $\mu A_{k} \left( {x_{j} } \right)$ that setting a mapping $\mu A_{k} \left( {x_{j} } \right):U \to \left[ {0,1} \right]$ characterizes $A_{k} \left( {x_{j} } \right)$ as the fuzzy set:

$$A_{k} \left( {x_{j} } \right) = \left\{ { \mu A_{kj} \left( {x_{j1} } \right)/x_{j1} , \mu A_{kj} \left( {x_{j2} } \right)/x_{j2} , \ldots , \mu A_{kj} \left( {x_{jb\left( j \right)} } \right)/x_{jb\left( j \right)} } \right\}$$

(39)

where $x_{j1} ,x_{j2} , \ldots ,x_{jb\left( j \right)}$ stand for the values that $x_{j}$ takes on.

The symbol $\mu_{X}$ will stand for the collection of membership functions describing the set of input variables $X$, that is,

$$\mu_{X} = \cup_{{\left( {k,j} \right)}} \left\{ {\mu A_{k} \left( {x_{j} } \right)} \right\},$$

(40)

Moreover, the pair $\left( {A_{X} ,\mu_{X} } \right)$ will stand for what is known as a fuzzy partition of the input domain $U$.

Respectively, the output variable $y_{s}$ with $s = 1,2,..,p$ associates to a set $L_{s}$, of linguistic terms $B_{k} \left( {y_{s} \left( {x_{1} , \ldots ,x_{n} } \right)} \right)$ namely

$$L_{s} = \left\{ {B_{k} \left( {y_{s} \left( {x_{1} , \ldots ,x_{n} } \right)} \right)|k = 1,2, \ldots ,c\left( s \right)} \right\}$$

(41)

Similarly,

$$B_{Y} = \cup_{s} L_{s}$$

(42)

will stand for the collection of linguistic terms that characterize $Y$. Also, for $B_{k} \left( {y_{s} } \right)$ we associate a membership function $\mu B_{k} \left( {y_{s} \left( {x_{1} , \ldots ,x_{n} } \right)} \right)$ such that the mapping $\mu B_{k} \left( {y_{s} } \right): V \to \left[ {0,1} \right]$ establishes the fuzzy set

$$B_{k} \left( {y_{s} \left( {x_{1} , \ldots ,x_{n} } \right)} \right) = \left\{ { \mu B_{k} \left( {y_{s1} } \right)/y_{s1} , \mu B_{k} \left( {y_{s2} } \right)/y_{s2} , \ldots , \mu B_{k} \left( {y_{sc\left( s \right)} } \right)/y_{sc\left( s \right))} } \right\}$$

(43)

where $y_{s1} , \ldots ,y_{sc\left( s \right)}$ denote the values that $y_{s} \left( {x_{1} , \ldots ,x_{n} } \right)$ acquires. Concurrently, we have the collection of membership functions tied to $Y$,

$$\mu_{Y} = \cup_{{\left( {k,s} \right)}} \left\{ {\mu B_{k} \left( {y_{s} \left( {x_{1} , \ldots ,x_{n} } \right)} \right)} \right\}$$

(44)

and concomitantly we may also say that the pair $\left( {B_{Y} ,\mu_{Y} } \right)$ sets a fuzzy partition for the output domain $V.$

Now, for $n > 1,$ we will consider the Cartesian product $L_{X}$ of fuzzy sets that designate the set of input variables $X,$ namely

$$L_{X} = \prod_{1}^{n} L_{j} ,$$

(45)

and use the symbol $m$ to denote its cardinality, that is,

$$m = {\# }\left( {L_{X} } \right).$$

(46)

This leads to

$$m = \prod_{1}^{n} q\left( j \right).$$

(47)

Moreover, the $i{\text{th}}$ element of the Cartesian product $L_{X}$. with $i = 1,2, \ldots ,m$ associates univocally with a $n - tuple$ of linguistic terms $A_{k} \left( {x_{j} } \right).$ That is, we have an ordering

$$i \to \left[ {A_{{k\left( {1,i} \right)}} \left( {x_{1} } \right), \ldots ,A_{{k\left( {j,i} \right)}} \left( {x_{j} } \right), \ldots ,A_{{k\left( {n,i} \right)}} \left( {x_{n} } \right)} \right],$$

(48)

where in the above $n - tuple\left[ {A_{{k\left( {j,i} \right)}} \left( {x_{j} } \right)} \right]_{1}^{n}$ the index $k\left( {j,i} \right)$ takes a value out of the set $\left\{ {1,2, \ldots ,q\left( j \right)} \right\}$. Moreover, the elements of the Cartesian product $L_{X}$ have form of a conjunction $Q^{i} \left( {x_{1} , \ldots ,x_{n} } \right)$, namely

$$Q^{i} \left( {x_{1} , \ldots ,x_{n} } \right):\left[ {x_{1} \,is\,A_{{k\left( {1,i} \right)}} \left( {x_{1} } \right),x_{2} \;is\,A_{{k\left( {2,i} \right)}} \left( {x_{2} } \right), \ldots , x_{n} \;is\;A_{{k\left( {n,i} \right)}} \left( {x_{n} } \right)} \right].$$

(49)

And it is by no way redundant to write down the equivalence

$$L_{X} = \cup_{1}^{m} \left\{ {Q^{i} \left( {x_{1} , \ldots ,x_{n} } \right)} \right\}.$$

(50)

Meanwhile, the set of membership functions associated to the conjunction $Q^{i} \left( {x_{1} , \ldots ,x_{n} } \right)$ will be denoted by means of the symbol $\mu Q^{i} \left( {x_{1} , \ldots ,x_{n} } \right)$ and formally expressed through:

$$\mu Q^{i} \left( {x_{1} , \ldots ,x_{n} } \right) = \cup_{j} \left\{ {\mu A_{{k\left( {j,i} \right)}} \left( {x_{j} } \right)} \right\}.$$

(51)

Moreover, for each value of the index $i = 1,2 , \ldots ,m$ we can consider a relationship

$$Q^{i} \left( {x_{1} , \ldots ,x_{n} } \right) \to R^{i}$$

(52)

being $R^{i}$ a rule that associates an antecedent $Q^{i} \left( {x_{1} , \ldots ,x_{n} } \right)$ to a consequent on the output value $y_{i}$. In the general fuzzy inference system $R^{i}$ takes the form

$$R^{i} :\left\{ {\begin{array}{*{20}l} {if} \hfill & : \hfill & {[x_{1} \;is\,A_{{k\left( {1,i} \right)}} \left( {x_{1} } \right),x_{2} \;is\,A_{{k\left( {2,i} \right)}} \left( {x_{2} } \right), \ldots , x_{n} \;is\;A_{{k\left( {n,i} \right)}} \left( {x_{n} } \right)]} \hfill \\ {then} \hfill & : \hfill & {\left[ {y_{i} \left( {x_{1} , \ldots ,x_{n} } \right)\;is\;B_{k\left( i \right)} \left( {y_{i} \left( {x_{1} , \ldots ,x_{n} } \right)} \right)} \right]} \hfill \\ \end{array} } \right\}$$

(53)

with the set of tied membership functions $\mu Q^{i} \left( {x_{1} , \ldots ,x_{n} } \right)$ given by Eq. (51).

For the case $n = 1$, we have only one input variable $x_{1}$. Hence, according to Eq. (37) the input space can be characterized by a number $q\left( 1 \right) \ge 1$ of linguistic terms, i.e., in this case we could have

$$A_{X} = \left\{ {A_{k } \left( {x_{1} } \right) |k = 1,2, \ldots .,q\left( 1 \right)} \right\}$$

(54)

and correspondingly from Eq. (40) $\mu_{X}$ becomes

$$\mu_{X} = \left\{ {\mu A_{1 } \left( {x_{1} } \right),..,\mu A_{q\left( 1 \right)} \left( {x_{1} } \right)} \right\}$$

(55)

and by virtue of Eq. (47) for this case we have

$$m = q\left( 1 \right).$$

(56)

And regarding the arrangement of the elements of the Cartesian product established by Eq. (48), in this case $i = 1,2, \ldots ,m,$ it suffices to advance a correspondence $i \to A_{{k\left( {1,i} \right)}} \left( {x_{1} } \right)$, where $k\left( {1,i} \right) = i$. Therefore, we can consider antecedent conjunctions

$$Q^{i} \left( {x_{1} } \right):[x_{1}\,is\,A_{i} \left( {x_{1} } \right)],$$

(57)

sustaining inferential rules $R^{i} :$

$$R^{i} :\left\{ {\begin{array}{*{20}c} {if} & {\begin{array}{*{20}c} : & {x_{1 } \;is\;A_{i} \left( {x_{1} } \right) } \\ \end{array} } \\ {then} & {\begin{array}{*{20}c} : & {y_{i} \;is\;B_{i} \left( {x_{1} } \right)} \\ \end{array} } \\ \end{array} } \right\}.$$

(58)

We also have

$$\mu Q^{i} \left( {x_{1} } \right) = \left\{ {\mu A_{i} \left( {x_{i} } \right)} \right\}.$$

(59)

In summary, for case $n \ge 1$ as stated by Eqs. (35)–(59) we may consider a general fuzzy inference system F as a construct including an application F: $X \to Y$ characterized by sets of fuzzy partitions $\left( {A_{X} ,\mu_{X} } \right)$ and $\left( {B_{Y} ,\mu_{Y} } \right),$ the set inference rules $R = \mathop {\text{U}}\nolimits_{1}^{m} \left\{ {R^{i} } \right\}$ and a defuzzification operator $D$ that associates to the fuzzy set $\left[ {y_{i} \,is\,B_{k\left( i \right)} \left( {y_{i} } \right)} \right]$ in Eq. (53) or a crisp value $y_{i}$ in $V.$

In the Takagi-Sugeno-Kang fuzzy inference system representation, we consider decision rules $R^{i}$ having an antecedent $Q^{i} \left( {x_{1} , \ldots ,x_{n} } \right)$ of the form given by Eq. (49) but with a consequent taking a crisp functional form $y_{i} = f^{i} \left( {x_{1} ,x_{2} , \ldots ,x_{n} } \right)$. That is, for $n > 1,$ in a TSK system we may envision inference rules

$$R^{i} :\left\{ {\begin{array}{*{20}l} {if} \hfill & : \hfill & {[x_{1} \;is\;A_{{k\left( {1,i} \right)}} \left( {x_{1} } \right),x_{2} \;is\;A_{{k\left( {2,i} \right)}} \left( {x_{2} } \right), \ldots , x_{n} \;is\;A_{{k\left( {n,i} \right)}} \left( {x_{n} } \right)]} \hfill \\ {then} \hfill & : \hfill & {y_{i} = f^{i} \left( {x_{1} ,x_{2} , \ldots ,x_{n} } \right)} \hfill \\ \end{array} } \right\}$$

(60)

for $i = 1,2, \ldots ,m,$ and with $\mu Q^{i} \left( {x_{1} , \ldots ,x_{n} } \right)$ given by Eq. (51).

Particularly, in defining a TSK system the function $f^{i} \left( {x_{1} ,x_{2} , \ldots ,x_{n} } \right)$ could have a linear form, that is,

$$f^{i} \left( {x_{1} ,x_{2} , \ldots ,x_{n} } \right) = p_{1}^{i} x_{1} + \cdots + p_{n}^{i} x_{n} + p_{n + 1}^{i} ,$$

(61)

where the numbers $p_{j}^{i}$ stand for parameters that can be empirically identified from data and being at least one nonzero. Nevertheless, in the general Takagi-Sugeno-Kang settings $f^{i} \left( {x_{1} ,x_{2} , \ldots ,x_{n} } \right)$ can take a nonlinear form.

An important concept tied to a TSK fuzzy model is that of the firing strength $\vartheta^{i} \left( {x_{1} ,x_{2} , \ldots ,x_{n} } \right)$ of the antecedent of a rule $R^{i}$. For the case $n > 1$ where the rules involve in their antecedents conjunctions of the form $Q^{i} \left( {x_{1} , \ldots ,x_{n} } \right)$ (cf. Eq. (62)) the firing strength $\vartheta^{i} \left( {x_{1} ,x_{2} , \ldots ,x_{n} } \right)$ is obtained through the algebraic product of involved membership functions $\mu Q^{i} \left( {x_{1} , \ldots ,x_{n} } \right),$ formally

$$\vartheta^{i} \left( {x_{1} ,x_{2} , \ldots ,x_{n} } \right) = \prod_{j = 1}^{n} \mu A_{{k\left( {j,i} \right)}} \left( {x_{j} } \right),$$

(62)

and the normalized firing strength $\varphi^{i} \left( {x_{1} ,x_{2} , \ldots ,x_{n} } \right)$ is defined through

$$\varphi^{i} \left( {x_{1} ,x_{2} , \ldots ,x_{n} } \right) = \frac{{\vartheta^{i} \left( {x_{1} ,x_{2} , \ldots ,x_{n} } \right)}}{{\mathop \sum \nolimits_{i} \vartheta^{i} \left( {x_{1} ,x_{2} , \ldots ,x_{n} } \right)}}.$$

(63)

The final output $y$ of the Takagi-Sugeno-Kang inference system, is the weighted average of all rule outputs, computed as

$$y = \mathop \sum \limits_{i} \varphi^{i} \left( {x_{1} ,x_{2} , \ldots ,x_{n} } \right)f^{i} \left( {x_{1} ,x_{2} , \ldots ,x_{n} } \right),$$

(64)

and identified as the output variable of the modeled system. The steps of the derivation of the formulae for the case $n = 1$ of a TSK fuzzy inference system are alike those establishing Eqs. (54)–(59) of the general fuzzy inference system and are elucidated in the methods section.

Appendix 4: Parameter Estimation for the TSK Fuzzy Model

In doing the parameter identification tasks, we bear in mind that the specification of the structure and the estimation of parameters of the TSK fuzzy model interrelates in such a way that none of them can be independently achieved without considering the other. Moreover, in the first stage of the involved identification tasks we made use of Substractive Clustering [8, 9], an efficient procedure endowing interesting advantages because it does not depend on optimization methods. Subtractive Clustering mainly relies on a measure of density of each data in the input space $X$. The goal is finding regions in the input space with high data densities. The data with the highest number of neighbors is selected as the center for a group. The data of the selected group that fall within a pre-specified fuzzy radius are then withdrew, and the algorithm keeps on searching for new data with the largest number of neighbors. This process is recursively applied until all the data are examined. Once Subtractive Clustering is completed, the number of decision rules can be known as each group associates to one of them. Moreover, the parameters tied to the membership functions characterizing the fuzzy sets of the antecedents of the rules are also estimated at this SC stage, and these are plugged into a recursive minimum squared (RLS) routine to obtain estimates of the parameters of the linear functions of the consequents of the rules. Following Jang et al. [20], we now explain the steps involved in the formalization of the RLS algorithm.

In the general least-squares problem, the output of a model $y$ is given by a linearly parameterized expression known as a regression function, formally

$$\varvec{y =} \varvec{\gamma_{1}} h_{1} \varvec{\left( u \right) + \gamma_{2}} h_{2} \varvec{\left( u \right) + \cdots + \gamma_{n}} h_{n} \varvec{\left( u \right)}$$

(65)

where $\varvec{u} = \left[ {u_{1} , \ldots ,u_{p} } \right]^{T}$ is the model’s input values vector, $h_{1} \left( \varvec{u} \right), \ldots ,h_{n} \left( \varvec{u} \right)$ are known functions of $\varvec{u}$, and $\gamma_{1} , \ldots ,\gamma_{n}$ called regression coefficients are parameters to be estimated.

Drooping the time variable for the sake of simplifying notation, we have that the general output of the TKS model of Eq. (7) becomes

$$w_{TSK} \left( t \right) = \varphi^{1} \left( a \right)f^{1} \left( a \right) + \varphi^{2} \left( a \right)f^{2} \left( a \right)$$

(66)

where $\varphi^{i} \left( a \right)$ are defined by Eq. (16),

$$f^{1} \left( a \right) = p_{1}^{1} a + p_{2}^{1}$$

(67)

$$f^{2} \left( a \right) = p_{1}^{2} a + p_{2}^{2} .$$

(68)

Replacing Eqs. (67) and (68) into (66) we have

$$w_{TSK} \left( t \right) = \varphi^{1} \left( a \right)\left[ {p_{1}^{1} a + p_{2}^{1} } \right] + \varphi^{2} \left( a \right)\left[ {p_{1}^{2} a + p_{2}^{2} } \right].$$

(69)

Then, rearranging we obtain that the regression function tied to model $w_{TSK}$ expressed in form (65) becomes

$$w_{TSK} = \varphi^{1} \left( a \right) p_{1}^{1} a + \varphi^{1} \left( a \right)p_{2}^{1} + \varphi^{2} \left( a \right) p_{1}^{2} a + \varphi^{2} \left( a \right) p_{2}^{2} a,$$

(70)

where $\varvec{a} = \left[ { a } \right]^{T}$ stand for the model’s input values vector, and $p_{1}^{1}$, $p_{2}^{1}$, $p_{1}^{2}$ and $p_{2}^{2}$ are the unknown parameters.

In order to get estimates for the involved parameters, we have to take into account that in the present settings the target system to be modeled involves an input-output relationship $a \to w\left( a \right)$ being $a$ the descriptor variable leaf area and $w\left( a \right)$ standing for the leaf biomass response. Therefore, we have a training data set composed of data pairs $\left( {a_{i} :w_{i} } \right)$, for $i = 1, \ldots ,m$ representing replicates of the addressed input-output relationship. Therefore, in order to identify the unknown parameters $p_{1}^{1}$, $p_{2}^{1}$, $p_{1}^{2}$ and $p_{2}^{2}$, in (70) we must fill in for each data pair $\left( {a_{i} :w_{i} } \right)$, into Eq. (70) in order to obtain the set of $m$ linear equations:

$$\left\{ {\begin{array}{*{20}c} {\varphi^{1} \left( {a_{1} } \right)p_{1}^{1} a_{1} + \varphi^{1} \left( {a_{1} } \right)p_{2}^{1} + \varphi^{2} \left( {a_{1} } \right)p_{1}^{2} a_{1} + \varphi^{2} \left( {a_{1} } \right)p_{2}^{2} } & = & {w_{1} } \\ {\varphi^{1} \left( {a_{2} } \right)p_{1}^{1} a_{2} + \varphi^{1} \left( {a_{2} } \right)p_{2}^{1} + \varphi^{2} \left( {a_{2} } \right)p_{1}^{2} a_{2} + \varphi^{2} \left( {a_{2} } \right)p_{2}^{2} } & = & {w_{2} } \\ \vdots & \vdots & \vdots \\ {\varphi^{1} \left( {a_{m} } \right)p_{1}^{1} a_{m} + \varphi^{1} \left( {a_{m} } \right)p_{2}^{1} + \varphi^{2} \left( {a_{m} } \right)p_{1}^{2} a_{m} + \varphi^{2} \left( {a_{m} } \right)p_{2}^{2} } & = & {w_{m} } \\ \end{array} } \right\}$$

(71)

In seeking for a solution to the above system we notice that it can be equivalently written in the concise form

$$BP = w,$$

(72)

where $B$ is the $m \times n$ matrix,

$$B = \left\{ {\begin{array}{*{20}c} { \varphi^{1} \left( {a_{1} } \right)a_{1} \varphi^{1} \left( {a_{1} } \right) \varphi^{2} \left( {a_{1} } \right)a_{1} \varphi^{2} \left( {a_{1} } \right)} \\ { \varphi^{1} \left( {a_{2} } \right)a_{2} \varphi^{1} \left( {a_{2} } \right) \varphi^{2} \left( {a_{2} } \right)a_{2} \varphi^{2} \left( {a_{2} } \right)} \\ {\begin{array}{*{20}c} \vdots \\ { \varphi^{1} \left( {a_{m} } \right)a_{m} \varphi^{1} \left( {a_{m} } \right) \varphi^{2} \left( {a_{m} } \right)a_{m} \varphi^{2} \left( {a_{m} } \right)} \\ \end{array} } \\ \end{array} } \right\}$$

(73)

and $P$ the $n \times 1$ vector of unknown parameters,

$$P = \left\{ {\begin{array}{*{20}c} {p_{1}^{1} } \\ {p_{2}^{1} } \\ {\begin{array}{*{20}c} {p_{1}^{2} } \\ {p_{2}^{2} } \\ \end{array} } \\ \end{array} } \right\}$$

(74)

And being $w$ the $m \times 1$ output values vector:

$$w = \left\{ {\begin{array}{*{20}c} {w_{1} } \\ {w_{2} } \\ {\begin{array}{*{20}c} \vdots \\ {w_{m} } \\ \end{array} } \\ \end{array} } \right\}.$$

(75)

The i-th row of the data matrix $\left[ {B \vdots \varvec{w}} \right]$ is denoted by $\left[ {b_{i}^{T} ,w_{i} } \right]$ and formally represented by,

$$b_{i}^{T} = \left[ {\varphi^{1} \left( {a_{i} } \right)a_{i} \varphi^{1} \left( {a_{i} } \right)\varphi^{2} \left( {a_{i} } \right)a_{i} \varphi^{2} \left( {a_{i} } \right)} \right]$$

(76)

Then, Eq. (72) is modified to incorporate an error vector $e$ in order to account for random noise or modeling error, that is,

$$BP + e = w.$$

(77)

Since $e = w - BP$ then $e^{T} e = \left( {w - BP} \right)^{T} \left( {w - BP} \right)$, and ${\text{if we let }}E\left( P \right) = e^{T} e$ we will have

$$E\left( P \right) = \mathop \sum \limits_{i = 1}^{m} \left( {w_{i} - b_{i}^{T} P} \right).^{2}$$

(78)

We call $E\left( P \right)$ the sum of squared errors and we need to search for $\hat{P}$, a characterization of the vector $P$, which minimizes $E\left( P \right)$. Moreover, the vector $\hat{P}$ is called the least-squares estimator (LSE) of $P$ and since $E\left( P \right)$ is in quadratic form, $\hat{P}$ is unique. It turns out that $\hat{P}$ satisfies the normal equation

$$B^{T} B\hat{P} = B^{T} w.$$

(79)

Furthermore, $\hat{P}$ is given by

$$\hat{P} = \left( {B^{T} B} \right)^{ - 1} B^{T} w$$

(80)

A $k$-order least squares estimator $\hat{P}_{k}$ of $\hat{P}$ defined by means of the expression

$$\hat{P}_{k} = \left( {B^{T} B} \right)^{ - 1} B^{T} w,$$

(81)

is a characterization of $\hat{P}$ that associates to $k$ data pairs taken out of the training data set $\left( {a_{i} :w_{i} } \right)$. Once we have obtained $\hat{P}_{k}$ we can get the succeeding estimator $\hat{P}_{k + 1}$ with a minimum of effort, by using use the recursive least-squares estimator (RLSE) technique, a procedure where the $k{\text{th}}\left( {1 \le k \le m} \right)$ row of $[ {B \vdots \varvec{w}} ]$, denoted by $[ {b_{k}^{T} \vdots w_{k} } ]$, is recursively obtained. In what follows, we will explain the formulism sustaining the RLSE method.

A new data pair $\left( {b_{k + 1}^{T} ;w_{k + 1} } \right)$ becomes available as the $\left( {k + 1} \right)$th entry in the data set, producing the $\hat{P}_{k + 1}$ estimate,

$$\hat{P}_{k + 1} = \left( {\left[ {\begin{array}{*{20}c} B \\ {b_{k + 1}^{T} } \\ \end{array} } \right]^{T} \left[ {\begin{array}{*{20}c} B \\ {b_{k + 1}^{T} } \\ \end{array} } \right]} \right)^{ - 1} \left[ {\begin{array}{*{20}c} B \\ {b_{k + 1}^{T} } \\ \end{array} } \right]^{T} \left[ {\begin{array}{*{20}c} \varvec{w} \\ {w_{k + 1} } \\ \end{array} } \right].$$

(82)

Further, in order to simplify the notation, the pair $\left( {b_{k + 1}^{T} ;w_{k + 1} } \right)$ will be symbolized by ($b^{T} ;w)$ and we also introduce the $n \times n$ matrices $H_{k}$ and $H_{k + 1}$ defined by means of

$$H_{k} = \left( {B^{T} B} \right)^{ - 1} ,$$

(83)

and

$$H_{k + 1} = \left( {\left[ {\begin{array}{*{20}c} B \\ {b^{T} } \\ \end{array} } \right]^{T} \left[ {\begin{array}{*{20}c} B \\ {b^{T} } \\ \end{array} } \right]} \right)^{ - 1}$$

(84)

or equivalently

$$H_{k + 1} = \left( {B^{T} B + bb^{T} } \right)^{ - 1} .$$

Then $H_{k}$ and $H_{k + 1}$ are related through

$$H_{k + 1} = \left( {H_{k}^{ - 1} + bb^{T} } \right)^{ - 1} ,$$

(85)

Therefore, using $H_{k}$ from Eq. (83) and $H_{k + 1}$ from Eq. (85), we explain why Eqs. (81) and (82) can be equivalently written in the form

$$\hat{P}_{k} = H_{k} B^{T} w$$

(86)

and

$$\hat{P}_{k + 1} = H_{k + 1} \left( {B^{T} w + bw} \right)$$

(87)

From (86) we have $B^{T} w = H_{k}^{ - 1} \hat{P}_{k}$, then replacing this result in Eq. (87) we get

$$\hat{P}_{k + 1} = H_{k + 1} \left( {H_{k}^{ - 1} \hat{P}_{k} + bw} \right).$$

(88)

Now, from Eq. (85) we have $H_{k}^{ - 1} \hat{P} = \left( {H_{k + 1}^{ - 1} - bb^{T} } \right)\hat{P}_{k}$, so replacing this result in the above expression we get

$$\hat{P}_{k + 1} = H_{k + 1} \left[ {\left( {H_{k + 1}^{ - 1} - bb^{T} } \right)\hat{P}_{k} + bw} \right],$$

then simplifying yields

$$\hat{P}_{k + 1} = \hat{P}_{k} + H_{k + 1} b\left( {w - w^{T} \hat{P}_{k} } \right)$$

(89)

Thus $\hat{P}_{k + 1}$ can indeed be recursively specified in terms of the previous estimate $\hat{P}_{k}$ and the new data pairs $\left( {b^{T} ;w} \right)$. Moreover, the current estimate $\hat{P}_{k + 1}$ is expressed as the previous one $\hat{P}_{k}$ plus a correcting term based on the new data $\left( {b^{T} ;w} \right)$; this amending term can be interpreted as an adaptation gain vector $H_{k + 1}$ multiplied by a prediction error $\left( {w - b^{T} \hat{P}_{k} } \right)$ linked to the previous estimator $\hat{P}_{k} .$

Calculating $H_{k + 1}$ as given by Eq. (84) is computationally costly and requires the adaptation of a recursive formula. From Eq. (85), we have

$$H_{k + 1} = \left( {H_{k}^{ - 1} + bb^{T} } \right)^{ - 1}$$

Applying the matrix inversion formula of Lemma 5.6 in Jang et al. [11] with $A = H_{k}^{ - 1}$, $B = b$, and $C = b^{T}$, we obtain the following recursive formula for $H_{k + 1}$ in terms of $H_{k}$:

$$H_{k + 1} = H_{k} - H_{k} b\left( {I + b^{T} H_{k} b} \right)^{ - 1} b^{T} H_{k} ,$$

equivalently,

$$H_{k + 1} = H_{k} - \frac{{H_{k} bb^{T} H_{k} }}{{I + b^{T} H_{k} b}}$$

(90)

In summary, the recursive least-squares estimator for the problem of $AP + e = y$, where the $k{\text{th}}{\mkern 1mu} \left( {1 \le k \le m} \right)$ row of $\left[ {B \vdots w} \right]$, denoted by $\left[ {b_{k}^{T} \vdots w_{k} } \right]$, is sequentially obtained, can be calculated as follows:

$$\left\{ {\begin{array}{*{20}c} {\hat{P}_{k + 1} = \hat{P}_{k} + H_{k + 1} b_{k + 1} \left( {w_{k + 1} - b_{k + 1}^{T} \hat{P}_{k} } \right)} \\ {H_{k + 1} = H_{k} - \frac{{H_{k} b_{k + 1} b_{k + 1}^{T} H_{k} }}{{I + b_{k + 1}^{T} H_{k} b_{k + 1} }}} \\ \end{array} } \right.,$$

(91)

where in establishing Eq. (91) we have taken into account Eq. (90) and the fact that we had formerly established the convention that in order to shorten the presentation the pair $\left( {b_{k + 1}^{T} ;w_{k + 1} } \right)$ would be symbolized by the simpler expression ($b^{T} ;w$).