Nonlinear Multivariable System Identification: A Novel Method Integrating T-S Identification and Multidimensional Membership Functions

Citation: Al-Hadithi, B.M.; Comina, M.; Jiménez, A. Nonlinear Multivariable System Identification: A Novel Method Integrating T-S Identification and Multidimensional Membership Functions. Appl. Sci. 2024,14, 6332. https://doi.org/ 10.3390/app14146332 Academic Editor: Angelo Luongo Received: 28 May 2024 Revised: 21 June 2024 Accepted: 24 June 2024 Published: 20 July 2024 Copyright: © 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). applied sciences Article Nonlinear Multivariable System Identification: A Novel Method Integrating T-S Identification and Multidimensional Membership Functions Mayra Comina 1,2,* , Basil Mohammed Al-Hadithi 2,3 and Agustín Jiménez 2 1Department of Energy Sciences and Mechanics, Universidad de las Fuerzas Armadas-ESPE, Sangolqui 171103, Ecuador 2Intelligent Control Group, Universidad Politécnica de Madrid, Centre for Automation and Robotics UPM-CSIC, C/J. Gutierrez Abascal, 28006 Madrid, Spain; [email protected] (B.M.A.-H.); [email protected] (A.J.) 3Department of Electrical, Electronics, Control Engineering and Applied Physics, Higher Technical School of Industrial Design and Engineering, Universidad Politécnica de Madrid, C/Ronda de Valencia, 28012 Madrid, Spain *Correspondence: [email protected] Abstract: In this paper, a new multidimensional Takagi–Sugeno (T-S) identification technique is proposed for multivariable nonlinear systems. In this technique, multidimensional membership functions are designed using concepts from solid mechanics. The design of membership functions is carried out in multidimensional space, defining the principal axes from the eigenvectors of the inertia matrix, and it has the characteristic of dividing the space into regions with the same inertia. These regions are analyzed to define the center of gravity for each rule. Illustrative examples of multivariable nonlinear systems, such as a thermal mixing process and a binary distillation column, are selected to evaluate the effectiveness of the proposed method. The proposed method is compared with traditional T-S identification that uses one-dimensional membership functions and shows a reduction in the relative identification error and the algorithm execution time. Additionally, the proposed method prevents rules from being positioned outside the system’s range, thereby avoiding the generation of unnecessary rules. Keywords: Takagi–Sugeno model; nonlinear systems; multidimensional membership functions; multivariable identification 1. Introduction System identification is the theory or art of constructing mathematical models of dynamic systems based on observed inputs and outputs. It involves creating models for unknown systems, simulating real-world behavior in cases where there is limited prior knowledge of the system’s structure. In this context, the identification of nonlinear systems is considered a challenging problem because it involves two main stages: selecting the model structure with a certain number of parameters and selecting an algorithm that estimates these parameters [1]. Given the complexity of modeling nonlinear systems, techniques like Takagi–Sugeno (T-S) are employed [ 2 ], providing a powerful tool for the identification and modeling of fuzzy systems due to their versatility, interpretability, and adaptability to a wide range of applications. These advantages make them valuable in control engineering, decision making, and other areas where the modeling and analysis of complex systems are required. 1.1. Fuzzy T-S Identification and Modeling In the literature, various approaches have been proposed for T-S fuzzy model identification. For instance, ref. [ 3 ] presents two novel learning algorithms for online T-S fuzzy Appl. Sci. 2024,14, 6332. https://doi.org/10.3390/app14146332 https://www.mdpi.com/journal/applsci Appl. Sci. 2024,14, 6332 2 of 46 model identification based on the Unscented Kalman Filter (UKF) and the dual estimation concept. These algorithms provide a more applicable parameter identification method by using the unscented transform instead of linearization, resulting in a more accurate propagation of the mean and covariance in highly nonlinear systems. Another implementation of T-S fuzzy model identification is found in [ 4 ], where parameter estimation of the nonlinear system is achieved by minimizing a quadratic cost. The study demonstrates that the system exhibits stable behavior and good transient response by designing a generalized T-S identification method along with an optimal state controller and an observer in each fuzzy rule. Similarly, the work by [ 5 ] employs a parameter weighting approach to optimize both local and global approximation in a T-S fuzzy model identification method. This methodology is pertinent to our study because it demonstrates how parameter weighting can significantly enhance the dynamic response and robustness of the system, ensuring zero steady-state error even under disturbances and modeling errors. Advanced methodologies for the identification and optimal control of nonlinear systems using generalized Takagi–Sugeno (T-S) models are presented in [ 5 , 6 ]. Both studies focus on precise parameter estimation by minimizing a quadratic performance index and employing parameter weighting techniques. These methods aim to optimize both local and global model approximation, resulting in systems that exhibit robust dynamic responses and effective damping. The application of evolving Takagi–Sugeno fuzzy models for nonlinear system identification is discussed in [ 7 ], showcasing the model’s ability to decompose input space into fuzzy subspaces, enhancing identification accuracy and computational efficiency. In [ 8 ], three subspace state-space algorithms are tested to select the appropriate algorithm that faithfully reproduces the dynamics of the real system with minimal estimation error. By combining T-S identification with enhancements in state-space techniques, a robust controller is achieved. Further advancements in T-S identification are illustrated in [ 9 ], where a tool is proposed for constructing T-S systems based on fuzzy c-regression models (FCRMs) and fuzzy clustering Gath–Geva. This method focuses on data clustering, particularly useful when domain expert data contain many errors, to construct fuzzy systems from experimental data. A novel data-driven methodology for constructing Takagi–Sugeno fuzzy models from experimental data is presented in [ 10 ], demonstrating the model’s effectiveness in approximating unknown nonlinear systems with fewer linearizations and higher control precision. Additionally, ref. [ 11 ] proposes a modified Gath–Geva fuzzy clustering algorithm for T-S model identification, optimizing model construction and achieving higher accuracy for nonlinear system identification. The use of fuzzy clustering techniques for system identification based on T-S fuzzy inference introduces innovative identification techniques. For example, ref. [ 12 ] defines a unidimensional Gaussian function as the antecedent of the rule, which showed better performance in all experiments and is attractive for areas such as data prediction. Another instance is the study presented in [ 13 ], which offers an extensive overview of evolving fuzzy and neuro-fuzzy methodologies in clustering, regression, identification, and classification. The researchers explore how these approaches can adapt and learn in real-time settings, progressively enhancing their understanding from incoming data. Notable applications in adaptive control and identification systems are discussed, illustrating the potential of these techniques to boost accuracy and efficiency across various practical scenarios. A method for Takagi–Sugeno fuzzy modeling using clustering algorithms to identify premise parameters is presented in [ 14 ]. The approach enhances generalization and approximation in dense input regions, achieving better control performance. The identification of nonlinear multivariable systems is further explored in [ 15 ], where the T-S fuzzy method based on backpropagation with gradient descent is used to model a quadcopter’s dynamics. The method demonstrated high identification accuracy (>98%) Appl. Sci. 2024,14, 6332 3 of 46 and computational efficiency. Similarly, ref. [ 16 ] applies T-S fuzzy models to describe the infectious dynamics of HIV with high accuracy using Gaussian fuzzy membership functions. Leveraging type-2 fuzzy T-S (T2-ETS) systems, a new technique for identifying nonlinear systems is presented in [ 17 ], balancing model complexity and prediction accuracy effectively. Moreover, ref. [ 18 ] introduces a combined methodology for the identification and adjustment of a Q-function, concluding that the proposed initialization reduces convergence problems and enables parameter optimization. The review of T-S fuzzy models for predictive control applications in [ 19 ] highlights the precision and adaptability of these models in controlling nonlinear systems, such as heat exchangers. The authors of [ 20 ] developed a method for identifying fuzzy systems based on entropy, demonstrating higher precision and efficiency in modeling the nonlinear dynamics of small helicopters. Addressing the control of uncertain nonlinear systems, hybrid neuro-fuzzy systems are proposed to overcome the drawbacks of individual fuzzy logic and neural network approaches. References to these methods are found in works such as [ 21 – 23 ]. An efficient procedure for nonlinear system identification using Takagi–Sugeno neuro-fuzzy models is proposed in [24], achieving high accuracy and robustness. The study in [ 25 ] integrates fuzzy stochastic configuration networks with Takagi–Sugeno inference, significantly enhancing the inference capability for complex nonlinear systems. Additionally, ref. [ 26 ] introduces type-2 Takagi–Sugeno fuzzy neural networks, showing the superior handling of uncertainty and improved system modeling accuracy. 1.2. Multidimensional Membership Functions, Multivariable Modeling and Identification The design of membership functions is the most crucial step in T-S identification for multivariable nonlinear systems. In [4,5,27,28], membership functions are designed using histograms, but this method can place fuzzy rules in inappropriate locations, whereas multidimensional membership functions are designed using point mapping, thus reducing identification errors. Multidimensional membership functions are designed in [ 29 ] using normalized histograms. The transformation from probability to possibility is resolved using two methods: the bijective transformationmethod[ 30 ]andthe uncertainty conservationmethod [ 31 ]. This methodassumes that the degree of membership is the same as the frequency of occurrence. A bidimensional synthetic Gaussian distribution is used for the probability-to-possibility transformation. In [ 32 ], multilayer feedforward neural networks are used to generate membership functions from labeled training data. One disadvantage of this method is that the shape of the membership function is unpredictable in regions where there are no training data. Membership functions are defined by linearly interpolating the membership degrees of characteristic points [ 33 ]. To achieve efficient interpolation, Delaunay triangulation of characteristic points is used, and interpolation is performed using barycentric coordinates. Composite and non-composite multidimensional fuzzy sets are proposed in [ 34 ]. A composite multidimensional fuzzy set is typically the result of combining two one-dimensional membership functions through union or intersection. Optimization techniques are also introduced to fine-tune parameterized membership functions for improved performance. In [ 35 ], multidimensional membership functions are presented in T-S fuzzy models for the modeling and identification of multivariable nonlinear systems using a genetic algorithm. This allows parameter convergence within a reasonable computation time and reduces identification errors with a smaller number of fuzzy rules. In [ 36 ], the enhancement of fuzzy modeling by directly utilizing multidimensional fuzzy membership functions to model complex, nonlinear systems is addressed. The authors demonstrate that this approach reduces the decomposition errors typically introduced by conventional methods, thereby improving model accuracy and performance. They employ radial basis function (RBF) networks to model the membership functions, showing significant improvements in transparency and interpretability. Appl. Sci. 2024,14, 6332 4 of 46 A novel method for identifying nonlinear systems using kernel functions is presented in [ 37 ]. The approach combines hierarchical identification principles with stochastic gradient algorithms, leading to improved parameter estimation accuracy. The study demonstrates that the proposed kernel functions, based on Gaussian membership functions, effectively handle the fuzzification and defuzzification processes, enhancing the overall performance of fuzzy logic control systems. The use of multidimensional membership functions in the fuzzy modeling of nonlinear systems is investigated in [ 38 ]. The authors propose a novel approach that maintains high model accuracy while reducing computational complexity. The results demonstrate significant improvements in model performance, particularly in handling complex nonlinear dynamics. An enhanced control strategy for nonlinear systems utilizing multidimensional membership functions is proposed in [ 39 ]. The authors show that this method significantly improves control accuracy and robustness by providing a more precise representation of the system’s dynamics. The study highlights the practical applications of this approach in various engineering fields. The application of multidimensional fuzzy membership functions in the identification of complex systems is explored in [ 40 ]. It is shown that this approach provides a more accurate and reliable model by capturing the intricate relationships within the data. The findings suggest significant potential for improving system identification processes in various domains. The proposed method in this paper is closely connected to previous Takagi–Sugeno (T-S) identification techniques that use one-dimensional membership functions, as presented in [ 4 – 6 , 27 , 28 ]. However, our approach introduces a significant improvement by designing multidimensional membership functions based on solid mechanics concepts, specifically using moments of inertia to define the principal axes and subdivide the data space into homogeneous regions. This approach not only reduces identification errors but also enhances the precision and computational efficiency of the model. The development of this method arose from the need to address the limitations observed in previous methods, where fuzzy rules were frequently positioned outside the system’s areas of influence, generating unnecessary computational load. Our process began with identifying these limitations through an exhaustive literature review and preliminary experiments that highlighted the need for a more robust and efficient approach. Thus, several iterations of multidimensional membership functions were designed and tested, culminating in the formulation presented in this work, which uses point cloud data from the system to avoid inappropriate rule placement The rest of this article is organized as follows. In Section 2, fuzzy T-S monodimensional modeling with classic formulation is described. In Section 3, multidimensional fuzzy T-S modeling is explained. In Section 4, the design of multidimensional membership functions is presented. In Section 5, illustrative examples of a multivariable thermal mixing process and a binary distillation column are explained to demonstrate the advantages of the proposed method. The results are presented in Section 6. The limitations as well as the contributions of this investigation are analyzed in Section 7. Finally, the conclusions are presented in Section 8. 2. Fuzzy T-S Monodimensional Modeling with Classic Formulation The proposed method utilizes generalized T-S identification to identify the nonlinear system in [ 5 , 6 ]. Identification is based on input–output data. The method estimates the parameters of the nonlinear system by minimizing a quadratic performance index using a parameter weighting method [ 2 ]. It models nonlinear functions as a set of difference equations with IF–THEN rules for an n-th-order system. Appl. Sci. 2024,14, 6332 5 of 46 The nomenclature ui is used for the inputs and yi for the outputs, considering p inputs. The method is based on identifying functions of the following form: f:Rp−→ R y=f(u1,u2,··· ,up) For identification purposes, we will also have m measurable variables (z1 , z2 , ··· , zm) , which can be inputs, outputs, or intermediate variables. Each IF–THEN rule S(i1...im) , for an n -th order multivariable system, with p inputs and q outputs, can be written as follows for each output: S(i1...im): If z1(k)is Mi1 1and z2(k)is Mi2 2and . . . and zm(k)is Mim mthen: yj(k+1) = a(i1...im) i0+a(i1...im) i1yj(k−1) + a(i1...im) i2yj(k−2) +···+a(i1...im) iniyj(k−ni) + p ∑ j=1 (b(i1...im) ij1uj(k−1)+ +b(i1...im) ij2uj(k−2) + ···+b(i1...im) ijniuj(k−ni)) (1) where y1(k) , y2(k) , ···yq(k) are the system outputs, and u1(k) , u2(k) , ···up(k) are the system inputs. The parameters ai0 , ai1 , ···aini are related to the system outputs, and the parameters bij1 , bij2 , ···bijni are related to the system inputs. As for the notation, j is the fuzzy variable index and ijis the fuzzy rule index associated with the fuzzy variable. The fuzzy estimation of the output is as follows: ˆ yj(k+1) = r1 ∑ i1=1··· rm ∑ im=1 w(i1...im)(z(k))ha(i1...im) i0+ a(i1...im) i1yj(k−1) + ···+a(i1...im) iniy(k−ni) + p ∑ j=1 (b(i1...im) ij1uj(k−1) + ···+b(i1...im) ijniuj(k−ni))   , r1 ∑ i1=1··· rm ∑ im=1 w(i1...im)(z(k))! (2) where w(i1...im)(z(k))= m ∏ s=1 µis s(zs(k)) (3) and the fuzzy rules weights w(i1...im) , and µis s(zs(k)) is a membership function of the fuzzy set Mij j . Let ns be a set of input–output system samples [x1k,x2k,..., xnk,yk] . The parameters of the fuzzy system can be calculated by minimizing the following quadratic performance index: J= ns ∑ k=1 (yk−ˆ yk)2=∥Y−XP∥2(4) where Y is an output vector, X is the input–output T-S matrix, and P is the T-S parameters vector. If X is a matrix of complete rank, the parameters of the fuzzy system are obtained as follows: J=∥Y−XP∥2=(Y−XP)T(Y−XP)(5) Appl. Sci. 2024,14, 6332 6 of 46 ∇J=XT(Y−XP)=XTY−XTXP =0 (6) P=XTX−1XTY(7) The problem exists if the membership functions are defined as triangular ones overlapped by pairs. In this case, the matrix X is not of full rank and thus is not invertible. An effective approach for solving this problem with low computational effort, based on the well-known parameters weighting method, is presented in [ 6 ]. This method can also be used for parameter tuning of the T-S model from local parameters obtained through the identification of a system in an operating region or from any physical input/output data. It is supposed that in this case, a first estimation of the parameters is available, which is obtained by the well-known least squares method. P0=p0 0p0 1p0 2. . . p0 nT(8) This first approximation can be utilized as reference parameters for all the subsystems. Then, the parameters vector of the fuzzy model can be obtained minimizing the following: J= ns ∑ k=1 (yk−ˆ yk)2+γ2r1 ∑ i1=1··· rm ∑ im=1 n ∑ j=0p0 j−p(i1...im) j2(9) J=∥Y−XP∥2+γ2∥p0−p∥2(10) J=   Y γp0−X γIp    2 (11) J=∥Ya−Xap∥2(12) where p0=P0P0··· P0T | {z } r1·r2···rm (13) In this case, the factor γ represents the degree of confidence of the parameters initially estimated. In a similar way to the previous equation, different weight factors of γ(i1...im) j can be used to each one of the parameters p(i1...im) j depending on the reliability of the initial parameter p0 jin the specific rule. The fuzzy discrete system can be described by the following IF–THEN rules for an n-th order system as a state model, as shown in [28]. S(i1...im): If z1(k)is Mi1 1and z2(k)is Mi2 2and . . . and zm(k)is Mim mthen: x(k+1) = A(i1...im) o+A(i1...im)x(k) + B(i1...im)u(k) y(k) = C(i1...im) o(k) + C(i1...im)x(k) (14) The discrete model of the fuzzy multivariable nonlinear system is represented as shown below: x(k+1) = Ao(k) + A(k)·x(k) + B(k)·u(k) y(k) = Co(k) + C(k)·x(k)(15) Appl. Sci. 2024,14, 6332 7 of 46 where the matrices A0(k) , C0(k) , A(k) , and B(k) are the fuzzy blends of the matrices A(i1. . . im) 0,C(i1. . . im) 0,A(i1. . . im), and B(i1. . . im)for a specific point. The effect of fuzzy rules on controller performance is related to the identification of the nonlinear system. By making the rules more fuzzy and choosing the appropriate universe of discourse, a precise identification of the system will be obtained, which will closely approximate the dynamic behavior of the nonlinear system. This, in turn, leads to the design of a robust and well-damped fuzzy controller that will meet the system requirements. However, increasing the number of fuzzy rules increases computational cost, so it is necessary to maintain a reasonable number of fuzzy rules. The process for determining the membership functions, and thus the number of fuzzy rules and fuzzy subsystems, is carried out through an iterative trial and error process. Generalized identification is tested with a fuzzy system configuration, and the identification error value is checked: e(k) = y(k)−ˆ y(k)(16) where ˆ yk is the estimated output of the fuzzy system described in (2). Ultimately, the resulting fuzzy system configuration is a trade-off between reducing the T-S identification error and maintaining a reasonable number of fuzzy rules. Equation (17) will be used to calculate the relative modeling error: erm(yi) = s∑k(yi(k)−ˆ yi(k))2 ∑k(yi(k))2=∥y−ˆ y∥ ∥y∥(17) 3. Multidimensional Fuzzy T-S Modeling It is proposed to use multidimensional membership functions (MDMFs) designed in multidimensional space to obtain better results than those obtained by using the traditional method. The approach proposed in this work is to use multidimensional fuzzy sets Mi where the fuzzy variable is defined as z= [ z1z2··· zf]t z∈Rf And the multidimensional membership function is denoted as µi(z). The T-S fuzzy model of the system with nmrules becomes Si1: If z(k)is Mi, then: y1(k+1) = ai 10 +ai 11y1(k) + ···+ai 1n1y1(k−n1+1) + p ∑ j=1 (bi 1j1uj(k) + ··· +bi 1n1uj(k−n1+1)) (18) ˆ y(k+1) = nm ∑ i=1 αi(z(k))hai 10 +ai 11y1(k) + ···+ai 1n1y1(k−n1+1) + p ∑ j=1 (bi 1j1uj(k) + ··· +bi 1n1uj(k−n1+1))#(19) αi(z(k)) = µi(z(k)) ∑nm i=1µi(z(k)) (20) Appl. Sci. 2024,14, 6332 8 of 46 4. Design of Multidimensional Membership Functions To minimize identification error using one-dimensional membership functions (1DMFs), the use of multidimensional membership functions (MDMFs) designed using the point map is proposed. To cluster the data, a geometric method is proposed using inertia axes that divide the space into regions, and then functions are assigned at the respective centers of gravity of the regions. The method used is explained in Figure 1. Figure 1. Steps for the design of multidimensional membership functions. The detailed explanation of each step of the proposed method is as follows. 4.1. Calculation of the Center of Gravity of the Central Rule The procedure begins with the point cloud yi , understanding that the central rule should be at the center of gravity, i.e., the mean of the point cloud. r1=1 n·∑ i yi(21) For example, in Figure 2, the area of influence of a system is shown where the variable yicorresponds to the cloud of points. Appl. Sci. 2024,14, 6332 9 of 46 Figure 2. Example 1—point cloud. Another example of a different system is shown in Figure 3, in which it is specified that the variable yicorresponds to the cloud of points. Figure 3. Example 2—point cloud. Then, the coordinates of the variable r1 , which corresponds to the center of gravity of the cloud of points of the system, are calculated. Figures 4and 5show the location of the center of gravity of the cloud of points for each example. Appl. Sci. 2024,14, 6332 16 of 46 The main objective is to keep Qs (outlet flow) and Ts (outlet temperature) at the desired values, which are considered the same as that of the tanks Q and T . The system is supposed to have a hot water with temperature T1 and cold water with temperature T2 . Both of them are considered constants. The two inlet flow rates Q1≥ 0 and Q2≥ 0 act as manipulable variables by means of two motorized valves whose dynamic model can be described as follows: Q1(s) = K1 1+T1se−τ1sU1(s) Q2(s) = K2 1+T2se−τ2sU2(s) (34) where •Q1and Q2: feed flow rates (m3/s); •T1and T2: feed temperatures of the tank, (◦C); •K1and K2: static gain of the system; •τ1and τ2: system time delay; •U1and U2: system inputs; • The “s” in Equation (34) represents the complex variable in the Laplace transform. It is also supposed that there are two measurements. It is further assumed that the measurement of the temperature T is carried out without delay. However, because of the location and technology of the flowmeter, it is supposed that there is a delay between the output flow and its measurement: Qm(t) = Qs(t−τs)(35) where •Qm: measured outlet flow, (m3/s); •Qs: outlet flow, (m3/s); •τs: time delay between the output flow and its measurement. Another assumption is that there are thermal losses proportional to the temperature difference between the inside and outside temperatures, which act as a disturbance effect. Thus, the model becomes Q1+Q2−Qs=Adh dt T1Q1+T2Q2−TQs+Kp(T−Te) = AdTh dt Qs=α√h (36) where •Kp: proportionality constant of thermal losses; •Te: environmental temperature, (◦C); •A: cross-sectional area of the tank, (m2); •α : proportionality constant which depends on the coefficients of discharge and the cross-sectional area; •h: height of liquid in the tank, (m). 5.1.2. Binary Distillation Column A distillation column is shown in Figure 14 with a feed flow Ff (mol/s) and its composition zfin the center of the column. Appl. Sci. 2024,14, 6332 17 of 46 Figure 14. Schematic diagram of the distillation column. There are upper and lower products, xD y xB . The system has two heat exchangers: a kettle at the bottom part that generates a rising steam Vf and a condenser in the upper part that cools the steam and generates the upper product and the reflux Ln . Part of the liquid is sent to the reflux tank where there is a liquid retention of mass MD0 (kg) with a composition xD. The reflux is pumped to the top plate NTat a rate of R0and ejected at the rate of Df. The distillation column has a plate structure to optimize the heat transfer between the liquid and the steam by maximizing the contact surface between them. In the base, the heavy products are removed at the rate of Bf with a composition xB and retention MB0(kg). Boiling steam is generated in a reboiler at the rate of V0. • Equilibrium phases In this work, a binary system (two components) is chosen with a constant relative volatility throughout the column and trays without losses (100% efficient); that is, the steam that leaves the tray is in an equilibrium state with the liquid in the tray. Yn=αXn 1+Xn(α−1)(37) A simple relation vapor–liquid equilibrium of (37) can be used for each tray. where –Yn: steam composition on plate n; –Xn: liquid composition of plate n; –α: relative volatility. • Hydraulic balance For the hydraulic balance, the dump in Figure 15 is considered, which allows the liquid to overflow from each tray in the distillation column. Appl. Sci. 2024,14, 6332 18 of 46 Figure 15. Schematic diagram of the distillation column. The molar flowrate of the outgoing liquid will depend on the fluid mechanics of the tray; that is, it is related to the liquid trapped in the plate. Ln=2p2g 3c·ρ·Lw(Mn−Mplate ρAplate )3 2(38) where –Ln: measured outlet flow (m3/s); –c: hydraulic discharge coefficient; –ρ: liquid density (kg/m3); –Lw: weir length (m); –Mn: accumulation of liquid of the plate, n+1; –Mplate: mass of the plate (kg); –Aplate: cross-sectional area of the plate (m2). To relate the accumulation of liquid in the tray ( Mn ) with the flow of liquid leaving the tray, (Ln) is used, as seen in (38). • Global mass balance and mass analysis by component. In a distillation column, a separation process stage is carried out in which, being a closed system, its mass will remain constant according to the law of conservation of mass. Each plate is considered as a stage in which the molar or mass component of the vapor and liquid flow is necessary to perform a mass and component balance. According to the schematic diagram (see Figure 16), the global mass analysis and mass analysis by component corresponds to Figure 16. Diagram of the dynamic balance of mass contact stage “n”. Appl. Sci. 2024,14, 6332 19 of 46 • Global mass analysis d dtMn=Ln+1+Vn+1−Ln−Vn(39) According to (39), each plate will receive a molar flow of liquid Ln+1 from the immediate superior plate and a molar flow of steam Vn−1 from the lower plate, and they will generate a molar flow of liquid of outlet Lnand a molar flow of steam Vn. • Mass analysis by component d dtMnXn=Ln+1Xn+1+Vn−1Yn−1−LnXn−VnYn(40) In the mass analysis by component according to (40), the respective concentrations of each of the flows must also be considered. Thus, the plate “n” receives the molar flow of the liquid Ln+1 with a concentration Xn+1 of the upper plate and a steam flow Vn−1 with a concentration Yn−1 of the lower plate. This in turns generates a molar flow of outlet liquid Ln with concentration Xn and a molar flow of steam Vnwith a molar composition of steam Yn. In conclusion, to determine the equations that govern the behaviour of a distillation column, it will be necessary to perform two types of analysis. The first one corresponds to a hydraulic analysis of the landfill through which the liquid flows from a higher plate to a lower one, and the second type corresponds to a balance of global mass and bycomponents. Due to the complexity of the model, the following assumptions are taken into account: 1. The heat trapped in each plate can be negligible. 2. The molar vaporization heats of the two components, distillate flow Df and sediment flow Bf, are approximately the same. 3. Heat losses from the column to the outside are negligible. These three assumptions lead to the following: V=V1=V2=... =Vn(41) The energy balance around each plate is not necessary. 4. The relative volatility αof the components is kept constant across the column. 5. Each plate has an efficiency of 100%; that is, the outgoing steam of each plate is in equilibrium with the plate liquid. 6. All the inflows and outflows of the tower are in liquid phase. 7. The feeding is achieved in a single plate. 8. There is no heat loss. 9. Steam accumulation is not considered throughout the system [42]. 5.2. Parameters for the Simulations This section details the parameters used for the simulation of the systems. 5.2.1. Multivariable Thermal Mixing Process For the simulation of the system, let us suppose the following process values: A=1 m2,α=1, T1=90 ◦C, T2=10 ◦C, Kp=0.1, t1=5 s, t2=4 s, t3=2 s, t4=3 s, τ1=1 s, τ2=0.5 s, τs=1 s The inlet flow rates are limited to 0 ≤Q1≤ 6 and 0 ≤Q2≤ 6, although in normal system operation, these extremes are not reached. The central operating point is shown below: Appl. Sci. 2024,14, 6332 20 of 46 u10 =2.0438, u20 =1.9563, Te=15 ◦C, Qs0=4 m3/s, T0=50 ◦C The sampling time is supposed to be ts=0.5 s. where •t1: time delay between the inlet u1and the measured outlet flow Qm, (s); •t2: time delay between the input u2and the measured output flow Qm, (s); •t3: time delay between the inlet u1and the outlet temperature Ts, (s); •t4: time delay between the inlet u2and the outlet temperature Ts, (s); •u10: initial condition of input 1; •u20: initial condition of input 2; •Qs0: initial condition of the flow rate, (m3/s); •T0: initial condition of the temperature, (◦C). 5.2.2. Binary Distillation Column The parameters of the binary distillation column are shown in Table 1. This table shows that the column operates with a feed flow of 2.5 kg/s and an input composition of 0.55. Table 1. Parameters of the distillation column. Variable Value Unit Nplates 20 - Nf10 - MD0500 kg MB0500 kg R00.75 kg/s V02 kg/s Ff2.5 kg/s zf0.55 - α2 - The initial reflux is 0.75 kg/s, and the initial steam flow is 2 kg/s with a thermal equilibrium constant of 2. The initial mass in the reflux drum and in the base-reboiler assembly is 500 kg in both cases. The hydraulic parameters of the column are required for the calculation of the mass in each plate for which the data in Table 2are used. Table 2. Hydraulic parameters of the distillation column. Variable Value Unit c 0.85 - g 9.81 m/s2 Lw1.225 m θ1.75 m hplate 0.06 m ρ997 kg/m3 5.3. Classical T-S Identification The implementation of the classical T-S identification method applied to the systems is shown below. Appl. Sci. 2024,14, 6332 21 of 46 5.3.1. Multivariable Thermal Mixing Process The design of one-dimensional membership functions is carried out using the method based on the histograms of the flow rate (Qm) and temperature (Ts). Figure 17a shows the histogram for flow rate Qm , and Figure 17b shows the histogram for temperature Ts. (a) (b) Figure 17. Histogram of Qm and Ts . (a) Histogram for flow rate Qm . (b) Histogram for temperature Ts . The fuzzy variables Qm and Ts are considered based on the histograms. The 1DMFs are defined as three triangular functions overlapped by pairs, as shown in Figures 18 and 19. Figure 18. Membership functions for flow rate Qm. Figure 19. Membership functions for flow rate Ts. Appl. Sci. 2024,14, 6332 22 of 46 From the data, the nine fuzzy rules shown in Figure 20 are obtained using the fuzzy inference method with 1DMFs. Figure 20. Fuzzy inference 1DMF. In Figure 21, it can be observed that the fuzzy inference with 1DMFs provides fuzzy rules whose results are outside the system’s range. Using these data places the rules in inappropriate points, leading to unnecessary computational resource usage. Figure 21. Central points of fuzzy rules using 1DMF fuzzy inference of thermal mixing process. Appl. Sci. 2024,14, 6332 23 of 46 5.3.2. Binary Distillation Column The method based on the histograms is used to design a one-dimensional membership functions of the bottom product (see Figure 22a) and the distillate product (see Figure 22b). (a) (b) Figure 22. Histogram of xB and xD . (a) Histogram for the bottom product xB . (b) Histogram for the distillate product xD. The 1DMFs for the bottom product and the distillate product are defined as three triangular functions overlapped by pairs, as shown in Figures 23 and 24. Figure 23. Membership functions for the bottom product xB. Figure 24. Membership functions for the distillate product xD. Appl. Sci. 2024,14, 6332 24 of 46 As seen in Figure 25, of the nine centers of gravity obtained with the classical T-S identification, seven are outside the system’s area of influence. Figure 25. Central points of fuzzy rules using 1DMF fuzzy inference of binary distillation column xD . 5.4. Multidimensional T-S Identification Using Five Rules Five rules are considered to test the multidimensional identification. These rules are bounded by the central point of the point map and by the central points of the half-planes defined by the inertia axes. 5.4.1. Multivariable Thermal Mixing Process In Figure 26a, the assignment of the central point of the central rule is obtained, and the inertia axes are shown in Figure 26b. (a) (b) Figure 26. Centers of gravity and level curves. (a). Center of gravity of the central rule. (b) Principal axes of inertia. In Figures 27 and 28, the subdivision of the point map into four zones delimited by the inertia axes is carried out. Appl. Sci. 2024,14, 6332 25 of 46 Figure 27. Case (a): zones 1–2 of the subdivision of the data points map. Figure 28. Case (a): zones 3–4 of the subdivision of the data points map. Then, the central points of the four zones are obtained (see Figure 29). Figure 29. Subdivision of the data points map. Appl. Sci. 2024,14, 6332 32 of 46 (a) (b) (c) (d) Figure 42. Level curves of the cases (a–d) for the binary distillation column. Finally, the multidimensional membership functions for cases (a)–(d) are obtained, as shown in Figure 43. (a) (b) Figure 43. Cont. Appl. Sci. 2024,14, 6332 33 of 46 (c) (d) Figure 43. Multidimensional membership functions of the cases (a–d) for the binary distillation column. 5.5. Multidimensional T-S Identification Using 9 Rules Multidimensional identification with nine rules is performed, utilizing the center of gravity of the four rules obtained from considering the half-planes defined by the inertia axes and the four rules from the four regions formed by the intersection of the inertia axes (see Figures 29,31, and 41). 5.5.1. Multivariable Thermal Mixing Process (e) Using for Qmand Ts−→ri=ci The centers of gravity of the nine rules are obtained using (25), and the multidimensional membership functions of the nine rules are observed in Figure 44. (a) (b) Figure 44. Case (e) for the multivariable thermal mixing process. (a) Centers of gravity of the rules. (b) Multidimensional membership functions (f) Using for Qmand Ts−→ri=2ci−r1 For the calculation of the centers of gravity of each rule for the flow rate and temperature, (26) is used. The centers of gravity of the nine rules and the membership functions are obtained according to Figure 45. Appl. Sci. 2024,14, 6332 34 of 46 (a) (b) Figure 45. Case (f) for the multivariable thermal mixing process. (a) Centers of gravity of the rules. (b) Multidimensional membership functions. (g) Different design criteria for each rule In this test, several design criteria for each rule are considered, as shown in Table 5. Table 5. Design criteria for the rules. Rule Criterion Central rule ri=ci Rules 2 to 5 ri=2ci−r1 Rules 6 to 9 ri=ci Thus, in Figure 46, the new centers of gravity and the multidimensional membership functions of the nine rules are shown. (a) (b) Figure 46. Case (g) for the multivariable thermal mixing process. (a) Centers of gravity of the rules. (b) Multidimensional membership functions. 5.5.2. Binary Distillation Column For cases (e), (f), and (g), the design criteria from Table 6are used. Table 6. Design criteria of the binary distillation column for 9 rules. Case xBxDEquation Subdivision Point Map (e) ri=ci(25) Figure 40 (f) ri=2ci−r1(26) (g) Central rule −→ ri=ci(25) Rule 2 to 5 −→ ri=2ci−r1(26) Rule 6 to 9 −→ ri=ci(25) Appl. Sci. 2024,14, 6332 35 of 46 The centers of gravity of the nine rules are obtained, as shown in Figure 47. (e) (f) (g) Figure 47. Centers of gravity of the rules of cases (e)–(g) for the binary distillation column. Finally, the multidimensional membership functions of the nine rules are obtained (see Figure 48). (e) (f) (g) Figure 48. Multidimensional membership functions of the cases (e–g) for the binary distillation column. Appl. Sci. 2024,14, 6332 36 of 46 6. Results of Multidimensional Identification To evaluate the effectiveness of the proposed method, three criteria are used: • Relative identification error; • Algorithm execution time; • Location of the centers of gravity of the rules. 6.1. Relative Identification Errors. The results of the relative identification errors obtained when applying the T-S identification algorithm using monodimensional and multidimensional membership functions are shown in Table 7. As observed in the obtained results, there is a lower identification error with the design criteria used in case (g). Therefore, case (g) is considered the best result obtained in this first parameter analyzed. Table 7. Identification errors. Test QmError TsError xBError xDError Least squares method 0.227 0.68 0.1449 0.1584 Traditional T-S with 9 rules 0.0033 0.0065 0.0945 0.0893 Case (a) 0.0024 0.0065 0.0056 9.1762 ×10−4 Case (b) 0.0027 0.0067 0.0058 9.5558 ×10−4 Case (c) 0.0029 0.0069 0.0078 0.0012 Case (d) 0.0026 0.0065 0.0052 9.7751 ×10−4 Case (e) 0.0025 0.0060 0.0044 8.9471 ×10−4 Case (f) 0.0024 0.0057 0.0039 8.0794 ×10−4 Case (g) 0.0024 0.0058 0.0036 7.9984 ×10−4 6.2. Algorithm Execution Time For measuring the algorithm’s execution time, 20 tests were conducted, and the results are shown in Table 8. Table 8. Algorithm execution times. Multivariable Thermal Mixing Process Binary Distillation Column Test Traditional T-S [µs] Multidimensional T-S [µs] Traditional T-S [µs] Multidimensional T-S [µs] 1 0.526674 0.340448 0.683674 0.482748 2 0.481436 0.453689 0.636436 0.597189 3 0.516763 0.451986 0.667763 0.600886 4 0.513475 0.430454 0.668475 0.582454 5 0.527426 0.383575 0.676426 0.528275 6 0.478236 0.420382 0.626236 0.566182 7 0.525911 0.410073 0.672911 0.556173 8 0.481825 0.440309 0.632825 0.589209 9 0.574699 0.451043 0.723699 0.595543 10 0.480659 0.430415 0.632659 0.573815 11 0.552393 0.407551 0.710393 0.556251 12 0.560221 0.460458 0.711221 0.610058 13 0.526787 0.414247 0.6766587 0.558447 14 0.502047 0.433955 0.650847 0.579555 15 0.48025 0.426617 0.63525 0.570717 16 0.49444 0.453739 0.64744 0.598239 17 0.522078 0.418289 0.673078 0.564589 18 0.534277 0.422008 0.686277 0.570908 19 0.518161 0.423949 0.667161 0.572749 20 0.484501 0.462671 0.632501 0.611071 Appl. Sci. 2024,14, 6332 37 of 46 From the obtained results in this test, it can be observed that there is an average execution time for the traditional T-S identification algorithm applied to the multivariable thermal mixing process of 0.514 µ s, and for the multidimensional identification, there is an average time of 0.426 µ s. For the case of the binary distillation column, the average execution time of the traditional T-S identification algorithm is 0.665 µ s, and for the multidimensional identification, the average time is 0.5732 µ s. Therefore, it is considered that the execution of the identification algorithm with multidimensional membership functions takes less time than the execution of the identification algorithm with monodimensional membership functions. 6.3. Location of the Centers of Gravity of the Rules The next parameter to be analyzed is the location of the centers of gravity obtained by implementing the T-S identification with multidimensional membership functions and monodimensional membership functions. 6.3.1. Multivariable Thermal Mixing Process For the multivariable thermal mixing process, the comparison between the centroids of the rules using the traditional T-S method and the proposed method of case (g) is shown in Figure 49. Figure 49. Comparison of the centroids of each rule using monodimensional and multidimensional membership functions for the multivariable thermal mixing process. As shown in Figure 49, the black indicators represent the centers of gravity of the rules obtained with the traditional T-S method, and it can be seen that some points are outside the Appl. Sci. 2024,14, 6332 38 of 46 system’s point map. In contrast, the centers of gravity of each rule obtained with the proposed method (red indicators) are within the system’s point map, verifying that our approach avoids placing the centers of gravity of the rules outside the system’s actual area of influence. 6.3.2. Binary Distillation Column By analyzing Figure 50, similar to the multivariable thermal mixing process, in the binary distillation column, it can also be observed that the centers of gravity obtained with the proposed method are within the system’s point map. Figure 50. Comparison of the centroids of each rule using monodimensional and multidimensional membership functions for the binary distillation column. 6.4. Compliance Analysis of the 3 Criteria When analyzing the results of the centroid locations obtained in cases (d) and (f), it was evident that despite achieving a reduction in identification errors, the centroids obtained are outside the system’s data points map, making it not an ideal solution for the problem at hand. However, in the remaining cases, all three conditions have been met, which means a reduction in identification errors has been achieved, and the centroids obtained are within the system’s data points map. It is concluded that the best system identification result is obtained with the identification algorithm using multidimensional membership functions with nine rules considering the centroid of the four rules coming from delimiting the half-planes defined by the inertia axes and the four rules of the four regions defined by the inertia axes. With the proposed methodology, there is a reduction in the identification error of 10.76% for temperature, 27.27% for the feed flow, 96.19% for the bottom product and 99.10% for the distillate product. 7. Discussion The comparison of the results obtained with the proposed method and previous studies are analyzed below. The findings and their implications are presented. The limitations of Appl. Sci. 2024,14, 6332 39 of 46 the work are discussed. Theoretical and practical contributions are highlighted, and future research directions are proposed. 7.1. Contributions of the Proposed Method After analyzing the results obtained from the proposed method, a comparative analysis of the reviewed articles on nonlinear system identification methods using multidimensional membership functions (see Section 1.2) is conducted compared with the proposed method (see Table 9). This specifies the relationship of the conducted research with each of the reviewed articles. Additionally, the unique contributions and advancements provided by this research are presented. Table 9. Cited articles—contributions of the proposed method. Article References Relation to the Proposed Method Contribution of the Proposed Method [4,5,27,28] Implements T-S identification, uses unidimensional membership functions, generates rules that are outside the system’s area of influence. Uses multidimensional membership functions to improve the accuracy of the fuzzy rule placement. [35] Uses genetic algorithms for the design of membership functions. The method is computationally intensive. This method reduces computational overhead by ensuring that the rules and their centers of gravity are strategically placed within the data map, thereby avoiding unnecessary computations. [36] Designs multidimensional membership functions using radial basis function (RBF) networks, presenting high complexity in optimization and parameter selection. Uses a physically-based approach (moments of inertia), facilitating model understanding and adjustment by engineers. [37] Designs membership functions based on Gaussian distributions, presenting a hierarchical approach to improve parameter estimation accuracy. Our proposal is specifically designed for nonlinear multivariable systems, providing a robust formulation for these complex environments, whereas [37]’s approach focuses mainly on improving parameter accuracy in a more general context. [38] The design and tuning of RBF networks can be complex, especially for systems with highly dynamic and unpredictable behaviors. Adequate parameter selection for RBF networks is crucial and may require a laborious iterative process. Uses moments of inertia to design membership functions, which is a more straightforward and less resource-intensive process. This reduces design complexity and facilitates implementation. [39] Designs multidimensional membership functions using interpolation methods and optimization techniques to improve nonlinear system control. Although it enhances nonlinear system control, it may not be as effective in reducing identification errors due to reliance on interpolation techniques. Using moments of inertia to design membership functions provides a more effective reduction in identification errors by better capturing the system’s distribution and dynamics. [40] Uses advanced interpolation and optimization techniques to design multidimensional membership functions. These techniques help better capture complex nonlinear relationships in the data. Although advanced interpolation and optimization techniques can improve accuracy, they may not ensure that rules are placed in the most appropriate intervals of the data space. Uses moments of inertia to ensure the precise placement of fuzzy rules, significantly reducing identification errors and improving model accuracy. 7.2. Comparison of Computational Complexity The computational complexity of the proposed method is compared with the identification algorithms reviewed in Section 1. Appl. Sci. 2024,14, 6332 40 of 46 From the computational complexity analysis in Table 10, it can be seen that the proposed method has certain advantages over the works cited in the Introduction Section 1 . Some of these advantages stem from the use of solid mechanics concepts. Table 10. Analysis of computational complexity. Method/Reference Method Description Computational Complexity Key Advantages Disadvantages Current Proposal T-S identification with membership functions based on moments of inertia Low, thanks to the simplicity of the formula and the efficiency in computing distances High accuracy, lower computational cost, differentiable, and easy to implement in real time May be less flexible for certain types of complex nonlinear data [3] UKF and dual estimation for T-S identification High (use of UKF can increase complexity) Higher accuracy in highly nonlinear systems, suitable for complex dynamics Intensive use of matrix calculations, which can be computationally expensive [4]T-S identification with quadratic cost minimization Medium (use of quadratic optimization) Good transient response and system stability Requires adjustments and computation time for quadratic optimization [5]Optimization of approximation with parameter weighting Medium–high (optimization with weighting) Robust dynamic response and zero steady-state error Greater complexity in selecting weights [7]Evolving T-S models for identification Medium (efficiency in input space decomposition) Improved accuracy and computational efficiency May require fine tuning of parameters [8]State-space subspace in T-S identification High (subspaces can be complex) Combines T-S with subspace techniques for robust identification Intensive use of subspace algorithms can be costly [9]T-S system based on FCRM models and clustering Medium (use of clustering and FCRM) Easy construction of fuzzy systems with experimental data May be less efficient for large datasets [10]Data-driven methodology for T-S models Low (fewer linearizations and high precision) High accuracy and control with lower computational cost May require more parameter tuning time [11]Modified Gath–Geva clustering for T-S Medium (clustering optimization) Higher accuracy in nonlinear system identification Additional complexity in clustering optimization [12]Clustering techniques for constructing T-S rules Medium (interpolation and clustering) Improved prediction accuracy and data handling Complex interpolation and potential scalability issues [14]T-S models with clustering algorithms Medium (clustering optimization) Improved generalization and accuracy in dense regions Complexity in parameter optimization [15]Quadcopter identification with T-S and backpropagation Low (use of only four membership functions) High accuracy and efficiency with low computational cost May be less suitable for more complex systems [16]T-S identification for HIV dynamics Medium (use of Gaussian functions) High accuracy in describing HIV infection dynamics Less efficient for systems outside the HIV context [17]Identification with T2-ETS systems Medium (balance between complexity and precision) Balance between complexity and prediction accuracy May require additional complexity adjustments [18]Combined method for Q-function in T-S High (use of LMI and Q-function adjustments) Reduction in convergence issues and optimization improvement Intensive use of LMI techniques can increase complexity [24]Neuro-fuzzy T-S models for nonlinear systems High (neuro-fuzzy model can be complex) High accuracy and robustness without complete system information Complexity in integrating neural networks and fuzzy logic [25]Stochastic configuration networks fuzzy T-S High (improvements in inference capability) Significant improvement in inference capability and identification accuracy Added complexity in stochastic configuration [26]Type-2 neural networks T-S for identification High (advanced handling of uncertainty) Better handling of uncertainty and improved modeling accuracy Complexity in designing and tuning type-2 networks 7.3. Computational Complexity To demonstrate the advantages of the proposed method in terms of processing time and computational resources, the results are compared with the nonlinear system identification methods reviewed in Section 1. Table 11 provides a clear and detailed comparison of computational complexity and allows for evaluating the performance of different methods in terms of processing time and required resources. Appl. Sci. 2024,14, 6332 41 of 46 It is important to note that the data were obtained from the referenced articles and serve only as a reference to provide an initial idea of the processing time and computational resources used when implementing the algorithms. To accurately compare the values, it would be necessary to implement all these algorithms on the same system. However, this is not feasible in this case, as each article works with different systems. Table 11. Processing time and computational resources. Method/Reference Processing Time (µs) Computational Resources Required Current Proposal (Multidimensional T-S) 0.426 Lower, which is due to the simplicity of multidimensional membership functions UKF and Dual Estimation [3]0.65 High, which is due to intensive matrix calculations and nonlinear transformations Quadratic Cost Minimization [4]0.60 Medium, which is due to the need for quadratic optimization Parameter Weighting [5]0.55 Medium, since weighting requires additional but not excessive calculations Evolving T-S Models [7]0.58 Medium, as there is efficiency in input space decomposition but requires fine-tuning State Subspace [8]0.67 High, since subspace techniques are complex and computationally costly 7.4. Theoretical Contributions The following are the theoretical contributions of the proposed method: • Introduction of Multidimensional Membership Functions (MDMFs) – The proposed method introduces a novel approach to design multidimensional membership functions using the inertia matrix derived from the system’s point cloud data. This approach ensures that the membership functions are well positioned within the system’s operational range, reducing the likelihood of generating unnecessary rules and enhancing the accuracy of the model. The use of inertia axes to subdivide the space into regions of equal “weight” and assigning membership functions to the centers of gravity of these regions are significant theoretical advancements. • Reduction of Identification Errors – By employing multidimensional membership functions, the proposed method significantly reduces identification errors compared to traditional one-dimensional membership functions. This is achieved through a more accurate representation of the system’s behavior, capturing the complexities of multivariable nonlinear systems more effectively. The proposed method demonstrates a notable reduction in identification errors for various test cases, including the multivariable thermal mixing process and the binary distillation column. • Efficiency in Computational Complexity – The proposed method offers a computationally efficient solution for system identification. The design of the membership functions using moments of inertia and the subsequent optimization process ensures that the computational cost remains low, which is crucial for real-time applications. The comparative analysis shows that the multidimensional identification algorithm has a shorter execution time than the traditional T-S identification algorithm, making it suitable for online control applications. • Robustness and Generalization – The multidimensional membership functions designed using the proposed method exhibit robust performance across different scenarios. The ability to generalize well in dense input regions enhances the control performance and robustness of the fuzzy models. This generalization capability is particularly beneficial for handling complex nonlinear dynamics in various engineering applications.