Model Based Predictive Control of Multivariable Hammerstein Processes with Fuzzy Logic Hypercube Interpolated Models

This paper introduces the Fuzzy Logic Hypercube Interpolator (FLHI) and demonstrates applications in control of multiple-input single-output (MISO) and multiple-input multiple-output (MIMO) processes with Hammerstein nonlinearities. FLHI consists of a Takagi-Sugeno fuzzy inference system where membership functions act as kernel functions of an interpolator. Conjunction of membership functions in an unitary hypercube space enables multivariable interpolation of N-dimensions. Membership functions act as interpolation kernels, such that choice of membership functions determines interpolation characteristics, allowing FLHI to behave as a nearest-neighbor, linear, cubic, spline or Lanczos interpolator, to name a few. The proposed interpolator is presented as a solution to the modeling problem of static nonlinearities since it is capable of modeling both a function and its inverse function. Three study cases from literature are presented, a single-input single-output (SISO) system, a MISO and a MIMO system. Good results are obtained regarding performance metrics such as set-point tracking, control variation and robustness. Results demonstrate applicability of the proposed method in modeling Hammerstein nonlinearities and their inverse functions for implementation of an output compensator with Model Based Predictive Control (MBPC), in particular Dynamic Matrix Control (DMC).


Introduction
The origins of Model Based Predictive Control can be traced to Model Algorithmic Control (MAC) [1] and Dynamic Matrix Control [2][3][4]. In this type of controller a model is used to predict a process' behavior, control actions are calculated aiming to minimize a cost function which is generally the quadratic error between a desired future setpoint and process' output. MBPC is advantageous in relation to other control techniques such as PID (Proportional Integral Derivative) and LQR (Linear Quadratic Regulator) controllers since it can consider actuator constraints, process' constraints and handle non-minimal and unstable processes as well as multivariable systems [5,6]. MBPC using Hammerstein models is still being investigated and applied both in numerical and experimental case studies, as depicted in Fig 1 (Query used in Scopus database: TITLE-ABS-KEY("Model Predictive Control" hammerstein)), showing the number of related publications on journals and conferences over the last decades and demonstrating scientific and academic interest in the research and development of this area in process control.
Many processes, in particular distillation columns and chemical reactors, can be modeled by Hammerstein models [7], where a nonlinear and memoryless static function precedes linear dynamics. An appropriate model of the Hammerstein nonlinearity is necessary to achieve adequate process control. Different approaches can be found in literature regarding this problem. [8] employs an Artificial Neural Network (ANN) to model the inverse function of the Hammerstein nonlinearity in a self-tuning configuration, demonstrating good results regarding representation of Volterra series expansion of Hammerstein model. Similarly, [9][10][11][12] present results and methods regarding applications of ANNs in modeling and control of Hammerstein processes. [13] employs two independent DMC controllers, embedding nonlinear equations to deal with static nonlinearities. [14,15] investigate the application of an output compensator based on the nonlinearity's inverse function with DMC, in particular, [15] proposes a decision rule in the event of multiple solutions for the inverse nonlinearity. [16] follows similar approaches but uses fuzzy models identified by recursive least squares (RLS), the resulting fuzzy Hammerstein models are either a single-input single-output (SISO) or multiple-input single-output system, which can be analytically inverted to obtain the inverse nonlinearity. [17] models a solid oxide fuel cell's nonlinearities using a multivariable fuzzy system and defines, similarly to [16], the inverse to be a straightforward analytical procedure since the resulting system is either SISO or MISO. In this paper we present an alternative solution to the control problem of systems with Hammerstein non-linearities. This solution can be applied to both monovariable or multivariable systems. Our solution is divided in two major steps. First, a general interpolator based on Takagi-Sugeno fuzzy logic is theorized and developed, named Fuzzy Logic Hypercube Interpolator, or FLHI, motivated by the necessity to adequately model static nonlinearities and its inverses. Second, FLHI is applied to modeling static nonlinearities and its inverses. These inverse models are employed as an output compensator for the predictive controller, resulting in an pseudo-linear system and allowing conventional linear control theory to be applied [8] to SISO, MISO or MIMO problems. Our proposal is as depicted in Figs 2 and 3, where e is the error signal between desired set-points y r and process outputs y, w Ã is the control action from DMC, u is the control action modified by FLHI considering the nonlinear static gain from block NL and w is the output from the static nonlinearity. In ideal situations w Ã = w, however, modeling uncertainties account for differences between the two signals.
Similar approaches can be found in literature such as [16,17], however, the following innovations are present and differentiate our work: i) fuzzy inference is defined on fuzzy logic operations, allowing changes such as choice of conjunction operator (t-norm) or choice of membership function to have major impacts on final results; ii) membership functions are used as kernel functions, allowing the interpolator to behave as nearest-neighbor, linear or cubic interpolator, to name a few; iii) the model's point cloud input space is divided in convex  regions and each region is projected (as in, mapped) to an unitary hypercube space, where interpolation occurs, this standardizes the input space from the fuzzy inference perspective and facilitates both interpolation and obtaining the inverse function; iv) inverse function is achieved by the solution of a nonlinear optimization problem for a known region, since the resulting fuzzy model is highly nonlinear depending on choices of conjunction operator or membership functions, as such, its inverse is not a "straightforward analytical procedure" as it is for other approaches in literature; v) multivalued (or multiset) inverse functions are adequately handled, as in multiple solutions are obtained if they exist, which allows for greater flexibility on control actions; vi) MIMO systems can be handled by our approach.
Results are presented for three study cases, a SISO system [18], a MISO system with uncoupled nonlinearities [19] and a MIMO system with highly coupled nonlinearities [20], presenting good results regarding control objectives such as reference tracking and minimization of control variation. Robustness considerations are also presented for cases where a mathematical model of the non-linearity is available.
The rest of this paper is organized as follows. DMC algorithm is presented for both monovariable and multivariable cases, further considerations are given in case of constraints as well as usage of output compensator for DMC. In what follows, foundations of the proposed Fuzzy Logic Hypercube Interpolator are presented. Then an overview and summary of the control problem using FLHI is given. Followed, are given considerations on robustness and performance metrics. Results are then presented for a SISO system, a MIMO system considering uncoupling and coupling of inputs and a MIMO system. In conclusion, final remarks are presented about the paper, the proposed method and future work.

Dynamic Matrix Control
Dynamic matrix control is one of the first model based predictive controllers, developed by Cutler and Ramaker [2][3][4]. Its internal model, the step response, is easily obtainable which allowed it to enjoy wide acceptance and industrial application, in particular chemical and oil [21] but also others such as automotive, food and aerospace [22]. Other advantages which contributed to its popularity are: applicability to multivariable systems; consideration of process constraints on inputs or outputs; prevention of excessive control actions; predictive reference tracking and disturbance rejection; to name a few [22,23]. DMC's Finite Step Response (FSR) internal model limits applications of the controler to open loop stable processes, however, alternatives are presented in literature [23,24] for unstable processes.
In what follows the DMC algorithm is detailed according to [6,23,25], first for the SISO problem and then extended to the MIMO problem.
SISO DMC Design. DMC aims to reduce future tracking error and control action increments by minimization of the cost function: whereŷ is the predicted process output j steps ahead given by a process model, y r is the desired set-point, N y is the prediction horizon, N u is the control horizon and λ is the move suppression factor. Process output prediction is given by the finite step response model: f is the free response, dependent only on past variables: Eqs (2) and (3) can be combined and rewritten in matrix form as: where:ŷ In Eqs (4) and (5) S T is an unitary vector with dimensions N y × 1, G is the dynamic matrix with dimension N y × N u and H is a matrix with dimension N y × N − 1: ; ð6Þ ðg N y þ1 À g 1 Þ ðg N y þ2 À g 2 Þ Á Á Á ðg NþN y À 1 À g NÀ 1 Þ : ð7Þ The functional Eq (1) can be rewritten in matrix form as: optimization of the control law is given by the minimization of this quadratic cost function in terms of control action increments. This is achieved by differentiating J with respect to control action increments vector u and equating to zero, i.e. @J/@Δu = 0. The resulting control law is given by: In practice, Eq (9) results in N u control action increments, however, only Δu(t) is used at each instant t. In the next instant t + 1 a new control action is calculated, this is known as sliding horizon control. Hence, only the first line of the gain matrix K dmc is needed, which helps reduction of computational effort.
MIMO DMC Design. For MIMO processes the effect of each input variable to each output variable is described by its FSR. Eqs (1), (2) and (3) are affected and must be rewritten to account for these extra variables. This can be accomplished in matrix notation of Eq (4), which helps in obtaining a low verbosity solution.
Considering a system with m inputs and n outputs, Eq (5) is rewritten to: Eqs (6) and (7) are rewritten in terms of G ij and H ij , the SISO matrices, for the i-eth output and j-eth input, as: ; ð12Þ Finally, the control law from Eq (9) can be applied considering a change of the vectors involved in prediction error: Constraints. When constraints are considered the optimum solution is no longer the analytical solution of Eq (9). In this case, iterative methods for quadratic programming are necessary [26] and the control problem can be rewritten as: A and b can be chosen to reflect limits on system variables such as, for example, control magnitude, process output magnitude or control increments [26].

Fuzzy Logic Hypercube Interpolator
The main goal in control by output compensator approach is precise identification and modeling of a process' non-linear characteristics so conventional linear control theory may be applied [8]. This problem motivated the creation of a general interpolator which exhibits desirable characteristics such as: modeling both a function and its inverse function; multivariable inputs and outputs; flexibility regarding interpolation characteristics; and, high computational efficiency.
In this section FLHI is presented according to it's working algorithm, which is separated in three different parts. In the first part, a user provided point cloud is verified for consistency and defines an internal model which is used to feed posterior calculations. All pre-calculations occur in the first part, which acts as a setup for the interpolator. In the second part, function interpolation is defined by a Takagi-Sugeno fuzzy inference system in an unitary hypercube space. In the third part, inverse function interpolation is defined as a root finding problem in hypercube space in terms of an optimization problem.
Interpolant Setup. At this initial stage, the expected user input is a point cloud and the output are regions for interpolations and respective hypercubes, the main component of FLHI interpolant. The point cloud is a set , where N is the number of points in the point cloud, xi = (xi 1 , . . ., xi m ) is a set of input coordinates of size m and xo = (xo 1 , . . ., xo n ) is a set of output coordinates of size n such that the generating function of the point cloud is a mapping f: < m ! < n . In this context, hypercubes are interpolation regions where input coordinates xi are mapped to an unitary space in the range [0, 1]. The algorithm for interpolation is defined in Algorithm 1.

Algorithm 1: FLHI Interpolant Setup Algorithm
Input : a pointcloud P Output : a set of regions // apply lexicographical ordering to the pointcloud according to input dimensions. 1 P = lexicographical_order(P) // initialize as empty set. A main characteristic of the point cloud is the distance between points for each dimension. A regular grid is determined by points which are equidistant across dimensions, that is in other words, with predetermined and uniform distances across dimensions. A semi-regular grid is determined by points with predetermined and non-uniform distances across dimensions. An irregular grid, i.e. scattered data, lacks structure or order regarding relative location of points.
Conversion from problem coordinates to unitary hypercube coordinates can be realized for a region by considering the base point as the null coordinate (0, 0, . . .) and adjusting each dimensions in all neighbor points in the region as either 1 if the neighbor's dimension moves away from the base point or 0 if it remains unchanged. This can also be achieved by mapping all points in a region to hypercube coordinates by: pointðiÞ:xiðjÞ ¼ pointðiÞ:xiðjÞ À base point:xiðjÞ xi stepðjÞ ; ð16Þ where m is the amount of input dimensions and xi step is the step size of a dimension, note base_point = point(1) in a region. Current proposal focuses on a point cloud forming either a regular grid or semi-regular grid. Irregular grid, i.e. scattered data, remains as future work but some considerations are presented for such scenario in this paper. Two main challenges follow irregular data: i) tesselating necessary regions for interpolation; and ii) mapping an irregular region to a regular hypercube.
Tesselation of regions for surface reconstruction from scattered data is an open research topic [27] and further investigations are necessary to find or develop a suitable algorithm for application in FLHI. This is made further challenging by the fact current methods in literature focus on triangulations that require 3 m points for a region, while FLHI is based on quadrangulations with requirement of 2 m points.
Mapping of an irregular region to a regular hypercube is feasible by current algorithms applied in finite element methods, such as projective transform or bilinear transform [28] mappings.
Interpolation. FLHI Interpolation can be separated in three major procedures, summarized in Algorithm 2, Algorithm 3 and Algorithm 4.
In Algorithm 2, a search occurs to determine which region produced by FLHI setup Algorithm 1 delimits desired interpolation coordinates, in problem domain, xi. Once the region is determined, desired interpolation coordinates must be mapped to the hypercube established by the region.

Algorithm 2: FLHI Interpolate
Input : a set of regions, a set of input coordinates xi Output : a set of output coordinates xo // find the region which contains the desired input coordinates. 1 region = find_region_input_contains(xi) // convert the desired input coordinates to hypercube space.
With all coordinates now in hypercube domain, bounded by [0, 1], interpolation occurs in the hypercube as presented in Algorithm 2. The main concept of FLHI is that each of the 2 m boundary points of the hypercube contain information about local geometry. Conjunction of information from all boundary points allows inference regarding true function value at any arbitrary position inside the hypercube. Thus, each of these points contribute with moment regarding whole hypercube area. Influence of each boundary point moment regarding any arbitrary position in the hypercube is inversely proportional to the distance between the boundary point and this arbitrary position.
Consider Fig 4, where two boundary points p1 and p2 are represented in the same input dimension x with unitary distance between each other. Each boundary point exhibits maximum, logical unitary, information at it's own position, for it is a sampled function value, and this information's contribution diminishes the further away from the point. Information contribution may be represented by a membership function that exhibits a maximum at the point's location and diminishes the further away from the point. Furthermore, it is unnecessary to define two different membership functions on a single dimension for it can be determined one is the logical complement of the other, that is: Previous logic can be extended to any number of dimensions in a logical hypercube, such that a point has multiple local membership functions lμ, one for each dimension of which it is composed. Global membership for a boundary point can given by the logical conjunction of all local membership functions for that point. Thus, global membership μ of a boundary point is given by: where T is the triangular norm (t-norm). The applied norm can be any of Godel, Lukasiewicz, Hamacher, Product, etc. In this paper the applied t-norm for all cases was the product norm. Finally, for each output dimension, an interpolated value can be obtained by first moment of area deffuzification: Membership functions can be defined arbitrarily and different interpolators may be obtained by appropriate choice of membership function. In this paper the following membership functions are explored: nearest neighbor, linear, cubic, lanczos and spline.
Inverse Interpolation. Inverse interpolation in FLHI occurs as described in Algorithm 5 and begins by searching which regions may output the desired interpolation set xo. If the desired output set xo fits in maximum and minimum output coordinate boundaries for a region, by the intermediate value theorem this region may produce the desired output. This process can be computationally sped up if, in FLHI setup, maximum and minimum output coordinate boundaries are determined for each region as a priori knowledge.
It is important to note that choices of t-norm and membership functions that lead to well defined logical hypercube space, where global memberships are bound in [0, 1], limit this process to the evaluation of maximum and minimum values of xo for each boundary point. However, in ill-behaved logical hypercube spaces, particularly for parametric membership functions such as cubic, a search must be performed to determine maximum and minimum values for each region.
When a region is determined as being able to interpolate the desired output coordinate xo, a root search procedure is performed in terms of xi. This is defined as the minimization of sum of squared residual errors between interpolation in this region and the expected output coordinates:

Algorithm 5: FLHI Interpolate Inverse
Input : a set of regions, a set of desired output coordinates xo Output : a set of input coordinates xi // initialize as empty set 1 xi = ; // check each region to see if its maximum and minimum outputs contain xo when continuous t-norm and membership functions are applied, interpolated space is continuous and existance of xo is guaranteed by the intermediate value theorem.

for each region in regions do 3 if region_output_contains(region, xo) then
// solve a minimization problem on variable xi using objective function, the sum of squared residuals. 4 x = minimize(objective_function, region.hypercube, xo) // convert hypercube coordinates to original problem coordinates. 5 x = convert_to_problem_coordinates(x) // add to set of solutions. 6 xi.add(x) 7 end 8 end Multiple regions may contain the desired output coordinates. As such, inverse interpolation is a multivalued function and may return multiple sets of solutions.

Control Algorithm Summary
This is the set of steps which summarize the proposed control approach in practice: • Setup a DMC controller with desired parameters N y , N u , N, λ and a step model from the linear block of the Hammerstein model; • Setup a FLHI interpolant with nonlinearity data from the static nonlinearity block of the Hammerstein model; • At each control instant, obtain process output and calculate the necessary linear control action of w(t) considering constraints on upper and lower bounds of nonlinearity; • Apply FLHI inverse interpolation with desired membership function on w(t) to produce the desired control signal u(t); • In case of multiple solutions from inverse interpolation, choose the one which minimizes control variation Δu(t); • Control loop is repeated as necessary.

Algorithm 6: FLHI Objective Function
Input : a set of test input coordinates xi, a set of expected output coordinates exo, an hypercube Output : sum of squared residuals error between exo and the evaluation of FLHI at xi // interpolate hypercube in desired position. 1 xo = interpolate_hypercube(hypercube, xi) // output sum of squared residuals.

Considerations on robustness of output compensation-Multiplicative gain uncertainty
Cancellation between the static nonlinearity function f and inverse model f −1 of FLHI is given by: such that in ideal conditions w(t) = w Ã (t). However, in practice, accuracy of FLHI's inverse model is not perfect but rather an approximation of missing information from the point cloud. This model uncertainty can be represented in Fig 3 by re-arranging the blocks as in Fig 5, where Δ m is an input gain uncertainty: where ideally the input gain uncertainty Δ m = 1. Output Compensation Stability Theorem. A system controlled by output compensation is asymptotically stable if the following necessary and sufficient condition is met: where GM and PM are respectively the gain margin and phase margin of the open loop gain H = DMC Á L, and |Δ m | 1 is the H-infinity norm of the input gain model uncertainty. Proof. A feedback, closed loop, system is asymptotically stable if the Bode-Nyquist stability criterion is met: Considering GM = 1/|H(jw pc )|, where jw pc is the phase crossing frequency and H is the open loop gain of an arbitrary system, Eq (24) becomes: Considering the input gain uncertainty Δ m is itself a static nonlinear gain that does not depend on frequency response, its worst case is given by its H-infinity norm ||Δ m || 1 and H = DMC Á L: Finally, by arranging terms Eq (23) is obtained: Measuring worst case model error. A definition of stability with output compensation control is given by Eq (23) considering worst case model error ||Δ m || 1 as a robustness metric in relation to gain margin. Model absolute relative error (MARE) can be represented by: then, worst case model error becomes a maximization problem: where ub and lb are respectively upper and lower bounds of input space x i . In practice the true nonlinearity NL is unknown but it is either mathematically or computationally modeled. In cases where only a point cloud from a real data set is available, this approach can be useful to measure the trade-off between a simple and a more complex model.
Given the locality nature of FLHI, originated from regions of interpolations, local optimization techniques are neither capable or satisfactory in solving Eq (30). Global search methods are necessary such as [29].

Performance metrics
In this section a performance metric is proposed to evaluate the effectiveness of different membership functions in FLHI models and its effects on set point tracking and control action. A fair assessment can be realized when the performance metric mimics the cost function Eq (1) of the model based predictive controller.
Set point tracking is evaluated by Integral Squared Error (ISE) of all outputs and output references: Control efforts are measured by the Integral Squared Variation of Control (ISVC): The last performance metric aims to mimic DMC's cost function Eq (1) and its purpose is to provide overall assessment of results: A final remark of caution is presented in regards to analysis of results. All results include an ideal case where only the linear process is controlled, disregarding nonlinearities. This ideal linear case is included to provide an estimate of optimal set point tracking and control variation, however, ISVC and J metrics for ideal cases consider the linear control signal w(t) instead of nonlinear control signal u(t), which is inexistent in these scenarios. Therefore disparities can be observed regarding ISVC and J metrics of ideal cases in contrast to nonlinear cases since different control magnitudes are involved, due to effects of static nonlinearities.

Results
In this section results are presented for three case studies in order to demonstrate the proposed method. The first study case regards a SISO system described in [18,30] where the nonlinearity is a fourth order polynomial. The second is a MISO system described in [19], its input nonlinearities are described by third order uncoupled polynomials. The third is a MIMO system described in [20,31], this system exhibits highly coupled input nonlinearities with exponential and polynomial terms.
For all case studies FLHI is used to model the system's nonlinear portion and its inverse for application in output compensation, then, DMC is designed considering the model's linear dynamics. Control action and process output dynamics are presented for all study cases as well as comparisons on the effects of different membership functions on DMC's cost function of Eq (1) and its relating performance indices ISE and ISVC. Results include the ideal scenario, where nonlinearities are ignored and only the linear system is controlled, as well as different membership functions such as nearest neighbor, linear, cubic, lanczos and spline.

SISO Study Case
A distillation column is modeled as a Hammerstein system and presented in [18]. In the original work [18], two models are given, one with a third order polynomial for the input nonlinearity and another with a fourth order polynomial. Both models exhibit a first order linearity.
A typical application of output compensation control [8] in this scenario would focus on the third order polynomial model since this can be trivially inverted and guarantees at least one real solution, being of uneven order. Despite being more accurate, the fourth order polynomial presents a problem in practical applications since its analytical inversion could lead to imaginary roots and an absence of feasible solution.
Our proposed method in this paper does not suffer from the problem of imaginary roots since model inversion occurs in the problem's universe of discourse. As such, parity of model order is not a problem for our approach. From [18,30], the fourth order polynomial Hammerstein model is: wðtÞ ¼ 1:04uðtÞ À 14:11uðtÞ 2 À 16:72uðtÞ 3 þ 562:75uðtÞ 4 ; Control results for FLHI with linear membership functions are presented in Fig 7. A comparison of the effects of different membership functions considering DMC's cost function is given in Table 1. For this study case, FLHI with nearest neighbor membership function fails to achieve reference tracking for all set-points and exhibits control ringing on the first set-point.
Robustness considerations are given for this study case in order to compare the effects of different membership functions on robustness. Robustness indices for this study case are GM = 10.3221, PM = 72.253 and maximum sensitivity M s = 1.226. Worst case model error considering MARE metric for each membership function is presented in Table 2, demonstrating nearest neighbor as the worst model and linear as the best. Stability criteria of Eq (23) is well met for all membership functions.

MISO Study Case
A Hammerstein system is proposed and employed in [19] for benchmarking a model identification method. This system, given in Eq (35), presents two inputs with uncoupled third order  polynomial nonlinearities. The linear subsystem presents coupling between inputs [32,33] according to its Relative Gain Array (RGA) in Eq (36).
Uncoupled model and control. When knowledge of uncoupling of nonlinearities is available it can greatly reduce computational needs. As such, in this first instance, this problem is modeled considering this knowledge. Control block diagram of this first approach is presented in Fig 8. A representation of the nonlinearity in Eq (35) Table 3. Regarding computational load, FLHI required a point cloud of 41 points and unitary hypercube dimension.
Coupled model and control. When coupling of nonlinearities is present or unknown, FLHI can be employed to model multivariable relationships. In this second instance, FLHI is applied considering both input variables are coupled even though in practice they are uncoupled. Control block diagram of this approach is presented in Fig 12. A representation of the nonlinearity in Eq (35) is presented in Figs 13 and 14, considering coupling and a nonlinear multivariable model. Process output and control actions are depicted  Table 4. Regarding computational load, FLHI required a point cloud of 1681 points and a two dimensional hypercube, respective to the system's input dimensions.
MISO Remarks. Set point was tracked in all cases and no ringing or abnormal control actions were observed. High ISE and J indices are explained by the very large jumps between set points.
Results were as expected regarding identical process response and performance indices for both uncoupled and coupled scenarios. A more complex model in this study case does not bring any benefit since the same behavior of uncoupling is modeled in both instances. FLHI's   increase in computational load in the second scenario is expected since the necessarily larger point cloud is a combination of both inputs and the internal hypercube mimics input dimensions, adding to the model's complexity.

Discussion
In this paper a novel method for modeling nonlinearities has been presented and applied to the problem of Hammerstein control using output compensation. The fuzzy logic hypercube interpolator, or FLHI for short, builds a fuzzy model based on point cloud data and allows model inversion. Model inversion enables FLHI to be applied directly as an output compensator, transforming the nonlinear control problem in a pseudo linear problem. Output compensation control, like in [8,16], is not related or anywhere similar to linearizing control [34,35].
Results are presented for a SISO process, a MISO process with uncoupled and coupled cases, and a MIMO process. These results include the ideal scenario, where only the linear system is controlled, and practical scenarios where the nonlinear Hammerstein system is controlled. In practical cases FLHI is applied using different membership functions such as Nearest Neighbor, Linear, Cubic, Lanczos and Spline. These results indicate the applicability of FLHI in both modeling of Hammerstein nonlinearities and output compensation, from its model inversion.
FLHI is currently limited to regular or semi-regular grid point clouds and injective data. Multivalued, i.e. non-injective, data and irregular grids are not automatically supported by the current method. Multivalued data can be used with FLHI but it is necessary for it to be manually separated in injective sets. Extrapolation is currently unsupported but the method can be trivially extended for it.
Future work includes, but is not limited to: i) support for irregular grids, as a possibility based on kd-trees; ii) support for multivalued data, using branch cuts; iii) study and analysis  of multiplicative uncertainties, modeling errors, created by FLHI in its application as an output compensator; iv) investigation of FLHI in modeling unknown Hammerstein nonlinearities; v) investigation of the applicability of FLHI in other control situations where model inversion is necessary; vi) study of different membership functions, in particular parametric ones similar to cubic; vii) study of different fuzzy t-norms, logical conjunction between membership functions.