Fuzzy model predictive control for small-scale biomass combustion furnaces

This work presents a fuzzy model predictive controller for small-scale grate furnaces based on a newly derived biomass combustion model. Several local linear controllers are designed for a selected number of operating points utilizing a gap metric. The resulting local predictive controllers are merged with membership functions to form a global nonlinear fuzzy control structure. The presented framework intends to improve the transient and steady state operation by applying an optimal control strategy with state estimation and to cover the entire operating range of the furnace. The open loop results of the introduced combustion model are parameterized and cross-validated with measured data from a test furnace. In order to find suitable parameters for the grey-box model, a local sensitivity analysis is conducted to contribute to an efficient parameter estimation process. Closed loop simulation results of the fuzzy model predictive controller, a linear model predictive controller and a PI control algorithm are presented and compared. Based on the performance of the proposed fuzzy controller, its application, advantages and disadvantages are discussed. Additionally, the impact of the different controllers on the formation of carbon monoxide is investigated based on estimation models from literature. The simulation results show that the fuzzy model predictive controller performs best in the considered categories.

• A state-of-the-art parameter sensitivity analysis is conducted for the model.
• An implementation of the nu-gap metric for model partitioning is presented.
• Fuzzy model predictive control for grate combustion furnaces is introduced.
• Carbon monoxide estimation models are applied to different control algorithms.

Keywords:
Biomass Modeling Fuzzy Predictive control Gap metric Emission A B S T R A C T This work presents a fuzzy model predictive controller for small-scale grate furnaces based on a newly derived biomass combustion model. Several local linear controllers are designed for a selected number of operating points utilizing a gap metric. The resulting local predictive controllers are merged with membership functions to form a global nonlinear fuzzy control structure. The presented framework intends to improve the transient and steady state operation by applying an optimal control strategy with state estimation and to cover the entire operating range of the furnace. The open loop results of the introduced combustion model are parameterized and cross-validated with measured data from a test furnace. In order to find suitable parameters for the grey-box model, a local sensitivity analysis is conducted to contribute to an efficient parameter estimation process. Closed loop simulation results of the fuzzy model predictive controller, a linear model predictive controller and a PI control algorithm are presented and compared. Based on the performance of the proposed fuzzy controller, its application, advantages and disadvantages are discussed. Additionally, the impact of the different controllers on the formation of carbon monoxide is investigated based on estimation models from literature. The simulation results show that the fuzzy model predictive controller performs best in the considered categories.

Introduction
During the last decades, the combustion of biofuels has been extensively researched in order to enhance the benefits of this renewable source of energy and to further exploit possible alternatives to fossil fuels. Not only the research, but also the practical application of biofuels has increased. In the year 2005, slightly more than 13.2% of the total energy produced in the European Union was obtained from renewable sources [1]. More than a decade later, in 2017 this share has already reached 22.9% of which 54.9% were obtained from the combustion of wood and other solid biofuels [2]. Detailed corresponding values for some selected European countries are given in Fig. 1.
Alongside the increased use of this renewable energy source, the release of CO 2 and other emissions like carbon monoxide, nitrogen oxides and particulate matter (PM) have increased as well and are influencing our environment and health [3]. Whereas CO 2 is the inevitable combustion product, other emissions like CO or NO x should be minimized. As a result, preventing the release of emissions is becoming of increased interest for most applications of biofuels. For the combustion of biomass, a well adjusted control strategy is an attractive way to decrease the formation of pollutants, to meet emission regulations and to further increase the performance of existing and new furnaces alike. The installation of additional filter is another way of preventing the release of pollutants, but it cannot influence their formation or increase the plant's efficiency directly. In this work, a fuzzy model predictive controller (FMPC) is introduced for small-scale grate combustion furnaces together with the underlying process model. The basic idea is to exploit an optimal control strategy that can handle nonlinearities, constraints and whose optimization can be extended by   1. The left bar represents the percentage of renewable energy to total energy produced in a country or region in 2017 [2]. The right bar refers to the share of energy obtained from the combustion of solid biofuels with respect to renewable energies.
additional control requirements, like emission estimation models or formation data as presented in [4]. A lot of unwanted gaseous emissions and PM result from incomplete combustion [5]. This often occurs during the transition between operating points, whereas a stationary operation is usually well adjusted. It is therefore necessary to provide a control strategy that is able to consider the transient behavior as well. Standard industrial control algorithms, like PID based controllers, are usually able to provide satisfying steady state results, but can lack performance during transition. Additionally, their implementation becomes more complex if several inputs or actuator limitations have to be handled. Therefore, more advanced model based control strategies are being implemented for the domestic and industrial combustion of biomass. Applications for model predictive control for large industrial furnaces are given for example in [6] for a furnace with moving grate, in [7] for a combined heat and power plant or in [8] for a municipal waste combustion plant. A different model based approach utilizing input-output linearization is given e.g. in [9] for a small-scale biomass combustion furnace whose design is comparable to the one considered in this work.
The general application of model based control algorithms however suffers from the well known problem of finding an adequate model which is simple enough for the controller design and the available hardware, yet complex enough to describe all relevant characteristics of the system. Depending on the field of interest, every modeling approach serves a different purpose and is therefore often unfit for most other applications. Although a lot of different furnace models exist, process modeling, emission formation and other problems are similar in different fields, as presented in [10,11] for the combustion of coal or in [12] for the co-combustion of coal and biomass. Combustion modeling based on distributed approaches as presented in [13] for a stratified downdraft gasifier and in [14] for fixed bed combustion focus on describing the complex thermochemical processes in detail, but are not suitable for control design. Other approaches as presented in [15] or in [16] are less complex, but sensitive equations are either used to describe the reaction kinetics or the heat transfer. In [6] an approach similar to the one in [16] was presented using a linearized Arrhenius equation, yet the overall model is still very detailed due to the presented furnace size and its application. Less complex process models based on simplified reaction kinetics as suggested in [17] are presented in [7,9,18]. In [7] a plant with a very specific structure is considered, whereas in [9] numerous experimental coefficients are introduced for the model. The model presented in [18] is simple and almost linear, but not suitable for the description of the operating range of the investigated furnace. As a result, a new process model is introduced in this work, which aims to use a reduced number of coefficients that can be estimated by nonlinear optimization. To emphasize the estimability of the parameters, a state-of-the-art sensitivity analysis based on the Fisher information matrix (FIM) is conducted [19] and the obtained parameters for the grey-box model are presented.
Based on the introduced process model, a network of constrained model predictive controllers is proposed, which is convoluted to form a fuzzy model predictive controller for small-scale furnaces. Linguistic and rule based Mamdani fuzzy systems [20] have been successfully applied in similar fields of research, for example to minimize N 2 O emissions [21], to estimate model parameters in combustion furnaces [22] or for the predictive modeling of biomass pyrolysis [23]. This work is based on Takagi-Sugeno (TS) fuzzy systems [24] however, for which Gaussian membership functions are applied favorably [25]. A survey about closed loop stability of these systems is presented in [26] and in detail e.g. in [27]. The local models for the linear controller design are obtained by splitting the operating range of the nonlinear process model into a number of linearized submodels. For this purpose a gapmetric approach is proposed based on the work presented in [28] and implemented as suggested in [29]. The selected -gap metric itself is based on H loop-shaping and combined with Gaussian membership functions reduces the number of independent parameters of the membership functions to one. The underlying MPC beneficially solves the problem of actuator saturation by a constrained optimization.
The combination of these methods yields an advanced overall control-scheme based on a nonlinear biomass combustion model for an existing plant with an application of the -gap metric. Closed loop simulation results of the FMPC for the small-scale biomass combustion furnace with a nominal load of 100 kW are illustrated in comparison to a global model predictive controller (GMPC) and a classic PID control concept from literature [30]. The different control designs are examined with regard to their stationary and transient performance and the resulting CO emissions are evaluated based on models for the combustion of wooden pellets [31]. This paper presents an approach on how to implement a fuzzy model predictive controller for a small-scale biomass combustion furnace and provides design methods for every step on the way. A byproduct of of the described procedure is the GMPC, which for itself is already an improvement compared to conventional furnace control. This work aims to introduce a fully method based state-of-the-art control concept for non-industrial furnaces and to evaluate and critically discuss the performance of less complex control algorithms in comparison.
This work is structured as follows: First, the underlying nonlinear furnace model is introduced. Next, the operating points for the local linear models are determined and the basic fuzzy formulation is presented. Then the time discrete model predictive control algorithm with the fuzzy convolution and state estimation is introduced. Finally, the model is validated for the investigated furnace and the simulation results for the controller are compared and discussed.

Furnace model
The introduced model of the furnace is based on physical principles complemented with experimental design approaches where missing information has to be obtained. The following sections present the basic furnace structure with two combustion stages, the nonlinear differential equations used for the controller design and a local sensitivity analysis for parameter estimation.

Process description
To obtain a suitable model for the controller, it is sufficient for small furnaces to use mass and energy balances without distributed parameters. Furthermore, by replacing complex reaction kinetics and heat transfer phenomena with reasonable simplifications, an appropriate process model can be defined by introducing only a manageable number of additional parameters. A schematic of the furnace is given in Fig. 2 where the massflows and temperatures available for modeling are depicted.
The inputs to the system are the fuel mass flow m f of wet biomass, the primary air from below the grate m pa and the secondary air m sa1 streaming into the freeboard above the grate where the second stage of the combustion takes place. The additional air inlet m sa2 allows for a more direct regulation of the flue gas burnout and therefore the flue gas composition. The combustion residues m ash and the flue gas m fg are process dependent mass flows leaving the furnace.
Owing to its size, the investigated plant is designed for district or domestic heating. The thermal power released during the combustion is distributed to the consumers via a water circuit. Therefore, the difference between the supply water temperature T sup and the return water temperature T ret of the heating circuit is the primary measurement of the plant power. The oxygen concentration of the flue gas O 2 contains valuable information about the process for feedback and control adjustments and indicates the current trade-off between performance and emission formation. If the O 2 concentration is small, efficiency is close to the furnace maximum but the concentration of CO is typically slightly higher than the achievable minimum [31]. High O 2 concentrations signal, that sufficient oxygen is available for a complete combustion but often indicate the begin of performance losses.
The freeboard temperature T fb and the temperature of the cold flue gas exiting the furnace T ex are important for a proper control strategy of the furnace as well, providing feedback about the combustion, emission formation, wear, soot and condensation in the exhaust system. The freeboard temperature is the main indicator for the thermal power released in the furnace, which should not be confused with the power of the water circuit. While the thermal power mainly reflects the amount of fuel that is converted during combustion, the water circuit heavily depends on the return water temperature T ret and the water mass flow m w as well. These two variables can be additional inputs to the furnace and for plant control, but are kept constant at the investigated furnace. The supply water temperature T sup is therefore the mixed product of these parameters. The exhaust gas temperature T ex can be considered mostly as a monitoring state, which serves as a lower temperature threshold for the exhaust system.

Solid mass balance
The current fuel mass m b of biomass on the grate can be obtained from with the mass flow of dry fuel m f,dry (in kg/s) onto the grate, m thd the thermal decomposition of the fuel (in kg/s) and the grate speed dependent disposal of combustion remains m ash (in kg/s). The dry fuel feed is given by where w H2O represents the mass fraction of the water content of the wet fuel mass flow m f (in kg/s). If the grate speed is constant and slow, such that a complete burnout of the combustible components on the grate is guaranteed, the relation can be applied with a ash as the ash content of the dry fuel. A simplified description of the thermal decomposition on the grate is presented in [17] as where k 1 represents the combustion rate constant (in 1/kg), m pa the mass flow of primary air (in kg/s) and m pa0 the load correction of the primary air (in kg/s).

Oxygen concentration in flue gas
The oxygen concentration in the flue gas after the combustion can be obtained from [30] (6) where m a is the total amount of air flowing into the furnace (in kg/s) including excess air and L min is the stoichiometric amount of air necessary for complete combustion given ideal conditions. Because laboratory conditions are not guaranteed in the furnace, the experimentally determined factor k 2 is introduced to correct L min towards measurements. The relationship presented in Eq. (5) yields good stationary and transient results for the oxygen concentration, but very fast dynamic effects are not represented entirely. Therefore the additional state R thd for the rate of change of the thermal decomposition m thd is introduced as with the experimentally determined time constant k Rthd . Eq. (5) is then modified to The additional term of the thermal decomposition rate allows a description of the observed fast dynamics of the oxygen concentration at the furnace.

Freeboard gas temperature
The temperature of the hot gas is obtained by an energy balance of the freeboard subsystem depicted in Fig. 2 where m g is the gas mass in the freeboard (in kg), c p,g is the specific heat capacity of the hot gas (in J/kgK), Q in is the enthalpy transported into the system with the air and the fuel (in W), Q comb is the heat released due to the combustion (in W), Q gas is the enthalpy of the hot gas leaving the freeboard to the heat exchanger (in W) and Q rad,fb are the radiation losses of the freeboard (in W). In detail, these terms are where c p,a is the specific heat capacity of air (in J/kgK), c p,s is the specific heat capacity of the solid fuel (in J/kgK), T amb the ambient temperature (in K), H L is the lower calorific heating value of the fuel (in J/kg), m fg is the flue gas mass flow (in kg/s), k 3 is an experimentally determined constant for the simplified radiation term (in m 2 K) and is the Stefan-Boltzmann constant (in W/m 2 K 4 ). The Stefan-Boltzmann law for thermal radiation requires the temperature to be of 4th order, but better results for the investigated furnace have been achieved with the simplified relation presented in Eq. (13). The enthalpy flow of the ash mass leaving the furnace is not considered in Eq. (9) due to its negligible magnitude.

Supply water temperature
The supply water temperature T sup for the heating circuit can be obtained by an energy balance of the heat exchanger. The thermal power of the transported water Q w is then part of the algebraic equation where Q ex is the enthalpy flow of the exhaust gas leaving the furnace (in W) and Q rec is the radiation recuperation term (in W), expressing the fact that not all the radiation energy of the freeboard is lost to the heat exchanger and the water circuit. The exhaust enthalpy flow is given by where c p,ex is the specific heat capacity (in J/kgK) and T ex is the temperature of the flue gas (in K). The radiation recuperation Q rec is again obtained by the simplified relation under the condition that < k k 4 3 is true and therefore Q rec cannot be a power source. The simplified radiation terms are beneficial, if the emissivity or the radiation surfaces of the furnace are unknown. With Eq. (14) linking the gaseous part of the furnace to the water part, the supply temperature T sup can be obtained from with the averaged heat capacity c p,w of water (in J/kgK), the water mass m w inside the heat exchanger (in kg), the water mass flow m w (in kg/s) and the temperature T ret of the returning water (in K).

Exhaust gas temperature
Assuming that no condensation occurs, the temperature of the gas at the end of the heat exchanger can be expressed with the enthalpy flow of the gas combined with a model for the heat transferred in the exchanger. A simplified approach is presented in [18] as where c p,ex is the constant averaged heat capacity of the gas (in J/kgK), m ex is the gas mass in the heat exchanger (in kg) and k R is a combined refractory dependent coefficient (in J/kgK). In combination with m fg , a load dependent heat transfer is obtained. However, Eq. (18) has a strong interdependency with Eq. (17), which restricts either the greybox estimation of the supply water temperature or the exhaust gas temperature. Since T sup is the major variable for control and T ex is not a primary target, the simple approach is suggested, which is basically a first order lowpass filter based on the freeboard temperature T fb corrected by the experimentally determined constant k 5 .

Parameter sensitivity and estimation
The nonlinear furnace model is represented by Eqs. (1)-(9), (13)- (17). The performance of the introduced model is presented and discussed in Section 3.1, after the parameters are obtained. The foundation of the parameter estimation process is described in the following. The states x m based on the first order differential equations of the model are comprised of the vector and the nonlinear continuous system can be expressed by The grey-box model also depends on the vector of yet undefined design parameters for the given model which have to be estimated based on measurements. The estimation is conducted with a particle swarm optimization (PSO) to overcome the obstacle of the results being dependent on the initial values. Introducing a high number or an unfavorable choice of design parameters however can lead to issues in the performance of the estimation process. Especially the reproducibility of the results can suffer due to nonunique solutions. Therefore, a state-of-the-art parameter sensitivity analysis based on the Fisher information matrix (FIM) as presented and discussed in [19] and in [32] is conducted in order to investigate the impact, magnitude and necessity of the parameters. The sensitivity analysis allows a priori statements about the local estimability of the parameters given experiment or measured data. The FIM can aid in grey-box estimation to identify problematic parameters or parameter combinations. These parameters typically show high variances, which means their contribution to explain the measurements is small. The FIM is based on the derivatives of the n x states with respect to the n p parameters in . Measurements are typically taken at a discrete time step denoted by k, which yields is the Jacobian matrix of the local parameter sensitivities. The Fisher information can be obtained as is a weighting matrix. The aim is to maximize the information contained in ( ), which is usually expressed by the Cramér-Rao inequality as var ( ) ( ) 1 (23) expressing that maximizing ( ) means to minimize the variance of the parameters. Because not all modeled states can be measured, the Fisher information is adapted based on the discrete state space representation given in [33] with where a system augmentation is introduced that elegantly resolves the partial derivatives of the parameters due to the expression A x k ( ) ( , ) m m and the Fisher information then is obtained as ( ). Based on the available measurements for O T T , , 2 fb sup and T ex , a set of initial parameters 0 is obtained which yields acceptable modeling results and allows an assessment of the range of the parameters. This range serves as a boundary for the applied PSO. The Fisher information is conducted for 0 and the parameter variances are evaluated. The results of ( ) 0 and the parameter estimation are presented and discussed in Section 3.1.

Submodels and fuzzy system
In order to design individual linear predictive controllers for different operating points, a local linear model network (LLMN) is required. In this section, the submodels for the model network and their locations within the operating range of the furnace are determined by a gap metric. The submodels used for control design are blended with Gaussian membership functions to form the overall nonlinear fuzzy structure.

The -gap metric
Various approaches exist to find sets of submodels, for example distributing linearization points equidistantly in the partitioning space, positioning using experience or with methods like local linear model trees [34]. In this work a gap metric is suggested which approaches partitioning from a controller point of view. Gap metrics are usually used for loop shaping in robust control, but some of their properties can also be applied for model spacing [35]. Based on the -gap metric introduced by Vinnicombe [28], the stability margin b P,C of the H norm is defined as if P C [ , ] is stable and as where P 2 is a perturbed version of P 1 . It is important to note that b 0 1 and 0 1 P,C By specifying a value for the margin and reformulating the relation from inequality (27), limits on the perturbation for plant P 2 with respect to P 1 can be found for which still holds. If the threshold depicted by is exceeded, a new submodel or linearization point is added, resulting in an overall closed loop transfer function based system description. A stepwise application of this method is presented in [29]. The -gap metric is conducted for all input-output combinations in order to find the largest gaps of the system. Discretizing the investigated operating range into 30 transfer function models with the margin = 0.25 yields the gap map depicted in Fig. 3 for the supply water temperature depending on the fuel mass flow. On the z-axis are the obtained gap values of the partitioning variable T sup and on the xand y-axis is the input m f normalized over its range, where 1 is equal to the stationary fuel amount for the nominal load of 100 kW and 0 equals the lowest operating point considered (30% load). The number of submodels determined by the -gap metric as depicted in Fig. 3 is = M 3 and is further discussed in Section 3.1. The supply water temperature is selected as partitioning variable to describe the transition between the local linear models. Alternatively, the total combustion air m a can be considered as independent variable instead of m f . This yields a different gap map and slightly different submodel locations, but because the amount of stoichiometric combustion air mainly depends on the fuel mass flow and not the other way around, this approach is discarded. Some of the following properties of the metric contain valuable information about the system: • If the perturbed plant model P 2 becomes unstable, the gap P P ( , ) 1 2 will become close to or exactly 1. This property can be used to identify problematic closed loop behavior, which is not the case for the derived model.
• Since the stability margin is fixed and the submodels are placed in the centers of the subregions, all local controllers are expected to perform equally well within their range.
• The method reveals which subregions can be represented by the same submodel, if their relative gap is smaller than , even if the subregions are not adjacent.
The -gap metric can also be applied to obtain the position for a single model or controller respectively. Based on Fig. 3, the location of the GMPC would be at 55% nominal load, because this is the location with the smallest maximum gap to all other operating points. In anticipation of the results however, the linearization point for the GMPC has to be close to the design point of the furnace for comparable results, which is at 100% nominal load. Therefore only the gaps between 75% and 100% are considered, which yields 90% nominal load as linearization point for the GMPC.

Membership functions
For the combination of the local MPCs, smooth and continuously differentiable transition functions are introduced to blend the submodels. If the transition is conducted with Gaussian membership functions and their conjunction is realized with the product operator a TS fuzzy system [24] is obtained, which is advantageous for the stability of discrete fuzzy controllers as shown ind [27]. Other membership functions are not considered in this work, but comparably good simulation results have been achieved as well.
Different means to blend submodels and to obtain global models or structures are available in literature. One possibility for example is to blend the outputs, which is a common approach in neuronal networks. The proposed FMPC structure in this work operates with a blended input vector. The following is based on the descriptions given in [25]. The input vector u from the fuzzy controller to the plant is composed of the weighted sum of the inputs of the = M 3 local linear MPCs as where x ( ) i par is the validity function depending on the partitioning variables x par , which is in this case but not necessarily only one state of the model, and u n i u is the input vector to the i-th subsystem. A necessary property of the validity function is that holds, which also guarantees that the combination of the inputs is still subject to the input constraints of the MPC. The validity or weighting functions themselves are obtained from as the Gaussian membership functions, where c ip is the location of the centerpoints and ip is the spread or standard deviation of the Gaussian. These two variables represent tuning parameters for the transition functions, but since the locations of the centerpoints correspond to the linearization points determined by the -gap metric, c ip is not available as degree of freedom. Due to the properties of the gap metric, local linear controllers located at the obtained points are expected to perform equally well [28]. This means, that the spreads of the Gaussians have to be equal too in order to preserve this property. Therefore, the remaining degrees of freedom for the parametrization of the membership functions for the M submodels are reduced to 1. Since p = 1, the product operator can be omitted in Eq. (33) and similarly i,p is replaced by M .

Control design
The intended network of fuzzy model predictive controllers consists of a number of individual MPCs, which are merged by the introduced membership functions. Each of these individual controllers is based on the state space representation of the submodels selected by the gap metric. Before the overall control structure for the network is presented, the integrated MPCs and their fundamental equations are introduced in this section.

Linear MPC with constraints
Because the design of the local MPCs is equivalent for all subsystems, the following formulation applies to all controllers with respect to their linearization point and the index i for the numbering of the submodels is omitted. The states x m , the controlled outputs y and the controllable inputs u are respectively. The primary air m pa and the secondary air m sa1 are aggregated, since they are not controllable individually for the given furnace. However, their ratio is fixed by geometry and a distinct treatment in the modeling equations is possible. The grate speed is preset in a way that only the ash content of the fuel is disposed, see Eq. (3), and therefore is no input to the system. The discrete state space representation of the local linear model for time step k is given by To eliminate steady state offsets, a set of integrator states is embedded into the system. For this matter a difference operation on both sides of Eq. (34) is conducted as shown in [36], yielding with the difference operators The matrices A and B are now given by The output vector can then be written in the compact form where × Q n n c y y is the weighting matrix for control errors and × R n n c u u the weighting matrix for control inputs. The optimal solution for the control increment U is obtained by setting the derivation of J c with respect to U equal to zero, yielding the control law subject to the constraints: Finally, the input vector u k ( ) of a local MPC is obtained by applying the receding horizon principle to U as where I m is an identity matrix and o m are N 1 c zero matrices mapping only the first n u elements of U to u k ( ) and discarding all other entries. The hard constraints on the outputs can be replaced by a set of soft constraints, if needed.

Fuzzy MPC with state estimation
The FMPC is obtained based on the local MPCs and the membership functions. An illustration of the fuzzy control structure is given in Fig. 4. Considering the individual controller operating points u o,i , the blended input vector u k ( ) of the fuzzy predictive controller is with the individual inputs u k ( ) Because not all states can be measured, an extended Kalman filter (EKF) is applied. The structure of the EKF for the state estimates x k ( ) based on [37] are with the additive process and measurement noise w k ( ) and v k ( ) respectively. Because the selected outputs correspond directly to modeled states, the nonlinear output function h (·) is replaced by the constant output matrix C m . The Kalman gain is updated as with the predicted estimation-error covariance matrix P k ( ) p and the covariance matrix of the measurement noise R Kal . The matrix P k ( ) p is obtained from with the updated state matrix A k ( ) m , the estimation-error covariance matrix P k ( ) and the covariance matrix Q Kal of the process noise. The covariance matrix for the next time step is obtained from where I is an identity matrix of size P. The control increment U from Eq. (39) for each MPC is then obtained by solving after replacing X k ( ) with the estimated state vector X k ( ).

Model parametrization and validation
This section presents the results of the parameter estimation, the introduced nonlinear furnace model and the local linear model network merged with the Gaussian membership functions. The results are obtained by means of measurements of the furnace using wooden pellets as fuel.

Parameter sensitivity and estimation
The variances of the initial model parameters 0 derived in Section 2.1.7 based on the Fisher information allow an a priori assessment of the parameters subject to the nonlinear parameter estimation. These results indicate which parameters are difficult to estimate with the given experiments, should be estimated separately or be removed from the model at all (i.e. they create a nonempty nullspace for ). This has already been conducted for the presented model. The remaining relevant parameter variances of 0 are depicted in Fig. 5 in the upper part. The obtained variances are upscaled by a large factor and then downscaled logarithmically, because otherwise only the last three presented values could be distinguished from zero in the same figure. This distorts the visible results slightly, but also highlights the critical parameters. Although the a priori variance for parameter k 1 is comparably small, it has to be removed from the estimation. This is due to the fact, that k 1 is an equilibrium parameter and only measurements of the mass on the grate can determine it definitely. The parameter is therefore fixed to yield the expected stationary amount of fuel on the grate, which allows the other time constants of the model to properly describe the remaining dynamic behavior. The parameter k Rthd is also excluded from the estimation, because it only describes the fast transient dynamics of the oxygen concentration, which are underrepresented in the available  measurements. Based on expert knowledge and observations, k Rthd is estimated to be about ten times as fast as k O2 .
The particle swarm used for identification is initialized with random parameter values uniformly distributed within their limits. Over a hundred independent estimation sequences are conducted in order to allow for a statistical meaningful a posteriori analysis and to compare the parameters to the a priori assessment based on the Fisher information. In the bottom part of Fig. 5, the scaled parameter values resulting from the nonlinear optimization are presented in the boxplot, after removing a small number of outliers. As predicted, the parameters k O2 and m g show the highest variance, whereas the other parameters converge very well to distinct values. The a posteriori parameter variance is therefore consistent with the a priori prediction of the Fisher information for the selected parameter subset. Since measurements were already available for estimation, the FIM can be applied for future design of experiments (DoE) in order to further define the parameters excluded from the optimization or reduce their variance in the estimation. The parameter values obtained for the investigated furnace are given in the Appendix.

Model validation results
To illustrate the performance of the model introduced in Section 2.1 and the local linear model network, open loop simulations are compared to measurements with a sampling time of Ts = 10sec obtained from the combustion of wooden pellets. In Fig. 6 the performance of the nonlinear and the local linear fuzzy model are presented on the identification data set. The model validation on new data is presented subsequently in Fig. 7. The data in Fig. 6 show a number of concatenated measurements of the system taken on different occasions at the investigated furnace. Various input steps have been applied to the system and the entire operating range is represented by different setpoints as well. To properly validate the nonlinear model and the fuzzy-LLMN, the model performance has to be tested on another set of measurements not available for training, which is presented in Fig. 7. The training and validation results are summarized quantitatively in Table 1 utilizing the Root Mean Square Error (RMSE) as a measure. The simulation results show that the nonlinear model is able to represent the open loop behavior quite well, especially T sup . Deviations occur at high oxygen concentrations and low temperatures, which are setpoints at which the furnace should not be operated for a long time due to the increased emission formation and inefficiency. It had to be expected, that not all dynamics are described equally well by the model due to the simplifications introduced in Section 2.1. However, some of the deviations between the model and the measured data are compensated by the extended Kalman filter. Additionally, not all deviations are equally severe for the MPC which is reflected in the entries of the weighting matrices Q c and R c used in Eq. (48).
The additional state for the oxygen concentration introduced in Eqs. (7)-(8) contributes to the fast observable dynamics of the process. Oxygen peaks occur at t = 1.8h in Fig. 6 and in Fig. 7 at t = 3.3h for example. The local linear model network shows in general a performance similar to the nonlinear model. Again especially T sup is captured well. The LLMN is convoluted with the membership functions determined for the controller and with the linearized models resulting from the margin = 0.25 for Eq. (29) for the -gap metric. Based on T sup three submodels for the operating range are specified. Their locations have already been presented qualitatively in Fig. 3 in the gap map and are given quantitatively in Table 2. The spread of the membership   , because based on the applied gap metric all controllers are expected to perform equally well within their designated range. The spread M and the centerpoints c i shown in Table 2 fully describe the applied membership functions µ which are depicted in Fig. 8 based on the constrained range of the supply water temperature. Because the remaining spread is equal for all Gaussians, they intersect exactly at half the distance between their centerpoints. In contrast to linguistic or Mamdani [20] fuzzy systems, no further defuzzification is necessary for the presented system. The partitioning variable T sup is transformed directly into the validity functions by Eqs. (31)- (33), which are utilized in Eq. (30) to form the weighted sum u of the local inputs for the fuzzy controller.

Closed loop simulation results
The simulation results of the fuzzy MPC presented in this section are compared to the global MPC and a PID controller based on the currently implemented control strategy. Both MPC concepts utilize the extended Kalman filter as presented in Section 2.3.2 for state estimation.

Simulation setup
The prediction horizon is set to = N 180 p and the control horizon to = N 60 c , which means that given the sampling time of Ts = 10sec N p covers half an hour. The comparison of these controllers illustrates the advantages of the FMPC and the GMPC and allows conclusions on when a PID control strategy can be considered. Since the derivative part of the implemented PID at the investigated furnace is zero, a PI controller is considered instead. The results for a series of steps in terms of nominal load of  Table 2.

Numerical comparison
To evaluate the controllers, the RMSE is applied to quantitatively indicate the control performance, see Table 3. Although steady state performance can hardly be evaluated by such a measure due to the integral part of the controllers, the general performance can be illustrated. The applied PI control concept is shown in [30], where the control error of T sup determines the inputs for the fuel mass flow m f and coupled by a constant factor (=k L 2 min ) the air mass flow m psa as well. The oxygen concentration has an additional control loop with the secondary air valve m sa2 as degree of freedom and can therefore achieve independent setpoints for O 2 to reduce emissions and guarantee a proper burnout of the combustibles in the flue gas. Although the control errors of the FMPC and the GMPC for the conducted series of steps are similar, all three controllers are able to eventually eliminate the steady state offsets over time. Comparing the dynamic properties of the step  Fig. 8. Membership functions µ for the three submodels depending on the supply water temperature T sup . The submodels refer to Table 2 and the centerpoints or linearization points are marked with the red-dashed lines.

Table 3
Comparison of the RMSE of the control error over the simulation time of the investigated controllers for the series of steps presented in Fig. 9 Table 4 and Table 5. For a fair comparison, all weighting matrices of the GMPC and the FMPC are identical. The FMPC however enables individual tuning for each local model, thus providing a higher potential for optimal settings. The dynamic specifications show a more clear distinction of the controller performances, highlighting the FMPC in almost all specifications. The differences between the control concepts of the PI and the MPCs are especially visible in the rise time t r for which the prediction makes a decisive difference. Additionally the FMPC shows slightly lower values for the overshoot than the GMPC.

Comparison of CO emissions
Closed loop simulation results are utilized to investigate the formation of carbon monoxide emissions at the furnace. For this purpose, the fuzzy-regression model presented in [31] for the CO estimation based on wooden pellets in small-scale biomass furnaces is applied. The controlled and constrained output variables O 2 and T fb are the inputs required for this black-box model. A new series of steps through the furnace power range is selected for which all effects of interest in the formation of CO can be qualitatively observed. The second profile in terms of decreasing and increasing the furnace power load is = Power [100, 75, 50, 30, 50, 75, 100]% and does not consider process and measurement noise in order to better highlight the different CO concentrations. Fig. 10 shows the simulation results for the estimated amount of CO normalized with the simulation time resulting in a quantitative indication of the concentrations. The predictive controllers are gaining momentum in the beginning in order to control the oxygen concentration, which explains the increased ppm levels at about t = 60min. The slow dynamics of the PI however become visible in the CO concentration as well, which peaks at t = 530min close to the end of the profile. This peak is explained by an unfavorable combination of the low temperatures in the freeboard but an already reduced oxygen concentration of the much faster second control loop. The estimated peak is even higher than the average CO levels at 30% load, where a stable combustion on the grate can barely be maintained. Fig. 10 shows independent of the underlying emission model, that in steady state all controllers eventually again will yield the same CO levels, but during state transition huge differences surface. This highlights especially the FMPCs ability to be adjusted to different transient dynamics, not only for control performance but also for emission limited combustion control. These two requirements however typically cannot be reconciled optimally by the same tuning, as they are opposing goals of a Pareto front [38] for combustion. Up to a certain degree, the FMPC offers the security to comply with the emission restrictions due to the output constraints and the soft transitions between operating points, therefore exceeding the capabilities of a common PI controller.

Applications of the FMPC
The introduced FMPC framework allows the application of easy to obtain local linear models, which are merged to form a global nonlinear control algorithm. As the presented FMPC is introduced mainly to account for the nonlinearities of the combustion process, other control aspects can be considered as well. As presented in [39] for example, a fuzzy controller is applied to limit the formation of different emissions during combustion. The necessary methodology for such an application at the investigated furnace is presented in [31]. Assuming that one MPC is sufficient and tuned for a fuel with a water content of 20% for example, a second MPC can be designed and tuned for the combustion of a similar fuel at 40% water content or a different fuel in general. The application of the fuzzy controller therefore can introduce various degrees of freedom, as each FMPC can have a different underlying model, tuning or fuel. For some of these applications, however, the main task will rather be to optimally merge the individual models.

Conclusion
A new combustion model for a small-scale biomass furnace has been introduced in the modeling section and open loop simulations have been compared to and validated with measured data. The introduced structure shows, that the nonlinear furnace model and the derived LLMN are capable of describing the most important characteristics of the combustion process. This is especially highlighted by the state considered most important in this work, the supply water temperature T sup , for which both models showed a high consistency for the available data. The presented modeling results translate to an average fit of about 80%, which is sufficient for a control application with state estimation.
The introduced model was profoundly analyzed with the Fisher information matrix to predict and support the estimability of the greybox parameters. The conducted runs of the PSO for the nonlinear parameter optimization delivered results within the expected variance for the parameters considered by the FIM. The -gap metric was applied to determine the linearization points for the network of linear submodels, which was combined favorably with Gaussian membership functions. The resulting nonlinear fuzzy structure covers the entire operating range of the furnace and requires only one additional parameter for tuning, which is very reasonable.
MPCs beneficially consider input and output constraints in their optimization and combined with an extended Kalman Filter, as shown in the control design section, provide a suitable approach. The FMPC was able to increase the stationary and transient performance of the GMPC, even if not always by far. Both predictive control concepts achieve about two times better stationary performances for the RMSE Table 4 Comparison of the average rise time t r in minutes of the different control algorithms for the series of steps presented in Fig. 9 Table 5 Comparison of the average percent overshoot PO of the different control algorithms for the series of steps presented in Fig. 9  than the PI controller. This performance indicator is however relaxed over time due to the integral part of the controller. Transient performance measures showed, that the FMPC outperforms the GMPC, but it is the prediction that brings the major improvements in the rise time by a factor of about four compared to the PI controller. The CO estimation results show, that it is important to keep the combustion parameters in ranges which provide a complete oxidation of the gaseous elements. The knowledge of these areas can be easily integrated into the constraints of a MPC, which can be key to achieve emission limited control. The advantages of the predictive controllers are obviously relaxed by the additional effort and their complexity. Therefore, the following suggestions are derived for the application of predictive controller for small-scale furnaces: • If mostly stationary control performance, very limited calculation power or easily manageable algorithms are of major concern, a classic PI(D) approach should be considered.
• If transient control performance, constraint control or general emission formation has to be considered, at least a linear model predictive controller should be applied.
• If the process is highly nonlinear, optimizing the performance or minimizing emissions is the major focus, the FMPC poses as a suitable candidate for small-scale furnaces.
This work introduced a comprehensive method based framework for fuzzy model predictive controller design together with a nonlinear furnace model and an application of the -gap metric and the FIM. The results indicate, that for the implementation of a predictive controller a single MPC might already be sufficient, depending on the furnace design and additional requirements. The next steps are the implementation of the presented control algorithm to the investigated furnace and analyzing the performance in comparison to the implemented PI controller.
The applied parameters obtained for the furnace model from nonlinear constrained optimization are: of which Q c is additionally weighted with the expected averaged steady state outputs ȳ and R c is weighted with the according averaged steady state inputs ū. The constraints for the inputs u k ( ), control increments U and outputs y k ( ) applied to the predictive controllers are presented in Table 6. To improve the calculation efficiency, the upper limit on O 2 and the lower limit on T fb can be omitted without affecting the results.