Design of Sliding Mode Control for Complex Nonlinear Models: A Numerical Approach

In this work, a numerical approach of formulating Sliding Model Control (SMC) is proposed to deal with complex nonlinear models of certain physical systems. Analytical formulation of SMC for such models is not possible as SMCs require algebraic manipulation of system equations, which is not achievable when the system model has strong nonlinearities. For such models, it is proposed to solve the SMC design using a numerical method that deals with numerical values instead of variables or functions in the equations. Mostly, SMC control law is dependent on the bounds of variables that are present in the mathematical expression of time derivatives of the sliding variable. To compute these bounds, a numerical approach is carried out using uncertainty bounds of the system’s parameters. Also, a numerical approach is used to compute time derivatives of nonlinear functions that are required in the formulation of SMC. Various forms of SMCs are investigated in this respect, including the basic first-order SMC and a second-order SMC. The proposed framework is generalized as well, making it compatible with a wide range of nonlinear models whose algebraic manipulation is not possible. A prototype example of a nonlinear model of a boiler is used as a proof-of-concept. The simulation results are promising and prove the efficacy of the proposed approach.


I. INTRODUCTION
Sliding mode controllers (SMC) belong to the class of controllers, which are robust and immune to parametric uncertainties of the model. The design of SMC control constitutes two phases; design of a stable sliding manifold that characterizes the system's dynamics in reduced dimensional form; design of a discontinuous control law that switches across the manifold in its neighborhood while maintaining the global reachability of sliding manifold [1]. The basic objective of the discontinuous control is to direct velocities of error trajectories towards the sliding manifold. The manifold itself is designed according to certain criteria, mostly using the Lyapunov stability theory, which enforces the reachability of manifold. Such reachability is to be strictly satisfied in order to achieve asymptotic or terminal stability of global error space.
Various implementations show the compatibility of SMCs with real-world systems [2], [3]. This is because, for practical The associate editor coordinating the review of this manuscript and approving it for publication was Bing Li . systems, many parameters have to be measured, and the uncertainties in the measurements can be troublesome for modeling and control purposes. To accomplish a robust controlled response of the system, it is always desired to use robust controllers like SMCs, where SMCs can show a high degree of robustness in the presence of uncertainties in the model.
The basic idea of sliding mode control emerged in the 1950s. The SMC in the original form constituted some weaknesses. For instance, the control input exhibits chattering in the vicinity of the sliding manifold, and although the manifold is reached asymptotically, it can sometimes lead to longer convergence times. Many modifications of SMC were done later in an attempt to tackle these shortcomings. Some of the works developed new kinds of manifolds and reaching laws so as to improve regulation and tracking performance. Others used Higher Order Sliding Modes (HOSMs) in sliding planes, where higher derivatives of the sliding variable are used to constitute sliding dynamics [4]. Under this category of HOSMs, Second Order SMC (2-SMC) controllers are the most researched and widely used in various practical VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ applications [5]. Many works have been reported that focus specifically on designing of unique forms of 2-SMCs [6], 2-SMCs with time delay [7], 2-SMCs with adaptive fuzzy control [8], etc. Super Twisting Controllers (STCs) under this category are particularly famous due to two reasons. One is, they are computationally efficient as they don't require the calculation of the derivative of the sliding variable, which is a necessity in other 2-SMCs formulations. Other is, they can work with the systems whose relative degree is one with respect to the sliding variable. Many articles prove the applicability of STCs with various kinds of physical systems. They have been analyzed, modified, and tested in conjunction with Neural networks [9], PIDs [10], backstepping controllers [11], and many others. There are many other forms of 2-SMC besides STC. However, a careful review of these SMC contributions determines that all the SMC literature is focused only on new designs of SMCs while disregarding the issue of the design complexity. Like other model-based control schemes, SMC uses algebraic manipulation of modeling equations to derive the control law, where the design process becomes very much cumbersome if the system models involve too many nonlinearities, lengthy equations, etc. Through the years, the complexity of new SMC techniques has grown, making them a tough candidate for control of highly nonlinear and complex models. One form of SMC, known as integral SMC, is used to improve the tracking performance of SMC to enhance robustness. Although integral SMC increases the order of sliding phase dynamics, which becomes equal to the system's dimension, the trajectory is made to start exactly at the sliding manifold, thereby eliminating the reaching phase [12]. Another line of SMC is dedicated to modifying the reaching dynamics by using hyperbolic terms [13], non-switching law [14], power terms [15], and many others [16] in control formulation. The target of such works is to develop unique kinds of sliding variable dynamic equations to improve the dynamics of reaching and sliding phases. Some methods attempt to improve the performance using different forms of sliding manifolds, for example, nonlinear manifolds [17], time-varying manifolds [18], moving manifolds [19], fuzzy integral sliding manifolds [20], etc. In the same line, PID manifolds have also been designed to ensure improved convergence to origin [21] or enhanced robustness for time-varying uncertain systems [22]. However, such manifolds increase the complexity of the SMC design process. This is because the analytical formulation of SMC requires lengthy algebraic computations, which cannot be carried out easily for complex models, which is further aggravated by increasing the complexity of the manifolds.
Another branch of SMCs is the terminal sliding mode, which tackles the shortcoming of asymptotic convergence of regular SMCs, by bringing error trajectories to manifold in a finite and shorter time. The terminal SMCs have a shortcoming that the control law may face the problem of singularity, which causes unbound growth of control inputs. This singularity issue has been dealt with different approaches [23]- [25]. However, all nonlinear manifold SMCs and especially terminal SMCs have the same problem, i.e., they use complex formulations of manifolds, which does not solve the problem of the design process for complex nonlinear models.
Artificial Neural Networks (ANNs) are also used in the SMC design process. The Radial Basis Function Neural Networks (RBFNNs) have been extensively used in this regard. In [26], Terminal SMC was improved by using an adaptive RBNN to estimate unknown function in the model of Permanent Magnet Linear Synchronous Motor, and the adaptive law was used to estimate the upper bound of the model's uncertainty. In another work [27], the authors used RBF in the design process of terminal SMC. With the proposed approach, the switching gain was adaptively adjusted using RBFNN to make the SMC robust and eliminate the chattering problem. In [28], RBFNN was used to model the dynamics of the robot manipulator, and the model was used in conjunction with SMC to execute a controlled response. All such works use ANN-based SMCs for modeling the unknown/unmodeled dynamics of the system or the sliding variables. However, they rarely touch the issue of complexity, and hence the proposed approaches are only compatible with dynamics models of limited complexity.
Regardless of all forms of SMCs, the mathematical designing of SMCs is considered as the most challenging part of it. The basic mathematical tools to design SMC are well developed and mature, but they are tedious or impossible to apply for complex nonlinear models. The complication increases with the increase in the number of nonlinearities and length of modeling equations. Particularly, the bounds and time derivatives of various functions become very hard to compute, as the complexity in the model increases. Such complexity issue is encountered in many models such as of a boiler, octarotor, quadrotor with tilting rotors, etc. Considering the case of the quadrotor, for instance, SMC has been used extensively with a simplified model of the quadrotor, but hardly has it been designed for the quadrotor with tilting rotors [29], [30]. This is because the additional four degrees of freedom afforded by the tilting mechanism increases the complexity of the model to a great extent (see [31] for example). The same is the case with robotic manipulators of sophisticated design or the octarotors [32], where the accurate model is too complex to design SMC unless a simplified model is used.
In this work, we aim to target this shortcoming when the design of SMC is not possible for complicated nonlinear models. In particular, we choose the nonlinear model of the boiler as the basis of our design framework. We use Astrom and Bell's model of boiler [33], which in fact, is the impetus of this work since all the previous works have focused on using a linearized model or simplified model of boilers for control design. Many practical system models pose the same complexity problem in the SMC design, due to which the proposed framework is also generalized in this work. Our target model [33] was derived using fundamental laws of conservation of energy and momentum. It is due to the model complexity that the previously reported works have used either a linearized model or a simplified model of boiler for the SMC control [34]- [36] etc. For instance, in [34], the SMC design was carried out using a PID sliding manifold for a linearized model of the boiler. In [36], the linear extended state observer (LESO) was implemented with the notion of compiling the total disturbances of the system as one state. The estimated disturbance was used in the formulation of SMC. Again, a simple model was used to design SMC as the level dynamics take extremely complicated mathematical form, which cannot be dealt with model-based control like SMC.
As a usual practice, the control law of the SMC sliding variable is mathematically dependant on the bounds and time derivatives of the variables involved in the expression of the sliding variable and its derivatives. For complicated models like boilers, the problem of analytical computation of the time derivatives and the variable bounds is either very complex or unsolvable. In this work, we propose a solution to this problem which is summarized as follows: 1) A numerical approach is proposed to solve the SMC design problem for complex nonlinear models with the aim of making the SMC design easily realizable without performing complex analytical computations. 2) The SMC design process constitutes the determination of the bounds of variables and the time derivatives of variables that are involved in the expression of the sliding variable. It is proposed that all such bounds and derivatives can be simply computed numerically. A deviation factor that depends on the parametric uncertainties is incorporated in these bounds, making them globally applicable around an operating point. 3) With the numerical computation of bounds and time derivatives, the SMC equations acquire a regular form for which any of the existing SMC algorithms can be applied with great ease.
The proposed method is applied and tested with various forms of SMCs on the complex model of a boiler system. It is proven that using a stable sliding manifold, the SMC control can regulate the level and pressure dynamics of a boiler under harsh action of steam disturbance, hence verifying that the proposed control technique works well in a simulation environment and has a potential of performing effectively in a real-time scenario as well. The paper is organized as follows. Section II is used to highlight issues behind SMC formulation of complex models and present the motivation of proposed work. Section III explains the model of the boiler and the associated complexities of the model. In Section IV, the boiler model is used as a build-up for numerical SMC formulation. Section V deals with few 2-SMCs control techniques. In Section VI, we propose numerical algorithms to compute time derivatives and bounds of variables involved in various SMC formulations. In Section VII, the general formulation of the SMC numerical framework is presented. Simulations and conclusions are presented in section VIII and section IX, respectively.

II. PROBLEM FORMULATION
Many control techniques involve strict use of system models to carry out the design phase of the control problem. SMCs are such controllers that highly rely on algebraic manipulation of modeling equations to formulate the control. In the SMC design procedure, the main aim is to make sliding manifolds an invariant set of global error space. The approach of the Lyapunov theory is used for the design process, which requires the time derivatives and the bounds of various variables involved in the SMC equations. It is observed that the major complexity lies in the analytical computation of the time derivatives, the variables bounds, and the typical form of the modeling equation, which is comfortably handled using available tools of SMC formulations.
In the SMC design, the complexity of various analytical computations can be avoided if the SMC design is carried out using a numerical approach. For the time derivates in SMC equations, we can consider a numerical approach of differentiation using forward or backward difference method. Such numerical differentiation can convert complex analytical functions y = h(x) into a linearized form y = C(x, t)x, a form that can be used easily for the SMC design using the available SMC tools The other problem is the bounds of complicated functions that can be addressed in a similar way. The idea is to use a numerical approach here to convert the uncertainty bounds of system parameters into the bounds of SMC variables. Such conversion can be carried out by discretizing the uncertainty bounds and using each discrete point to compute the numerical value of concerning functions in the sliding variable equation. The bounds of the functions would correspond to the maximum and minimum of the numerical values of the concerning functions, which are otherwise impossible to derive or compute using the analytical approach.
We implement the above procedure for a boiler system of the model [33]. For completeness, we present the boiler model first, and then we describe how a typical first-order SMC (1-SMC) can be used to control the state of pressure using Utkin's approach of SMC formulation. The level dynamics, in contrast to the pressure, are a complicated nonlinear function of boiler variables so we develop a numerical approach in an algorithmic form to design the SMC for such a complicated case of nonlinear output. In a similar fashion, we formulate a numerical approach for designing 2-SMC for the pressure control of the boiler system. Particularly, we numerically derive the time derivatives and bounds of SMC variables, which are used to formulate the 2-SMC for the pressure control system.

III. BOILER WORKING PRINCIPLE AND MODEL
A boiler is a device that is used to convert water into steam for various industrial applications. It is considered to be dangerous equipment since it operates at extremely high pressures and temperatures, which can be life-threatening if not adequately controlled. VOLUME 8, 2020 Many complications are associated with controlling a boiler system. Among them, the most intricate challenge is model-based control. It is because the model of the boiler system has many nonlinearities; analytical nonlinearities (variable multiplications, logarithmic and high power terms, etc.) and the interpolation nonlinearities of the steam table. Very effective control is required to be implemented in order to keep the states of the boiler under safety level. The controller should be versatile enough to perform all the essential tasks of stabilizing the unstable dynamics, performing the setpoint tracking for all outputs, and controlling the outputs without affecting the MIMO couplings that exist between all the inputs and the outputs.
In this work, we use Astrom and Bell [33] dynamic model of boiler, which uses four states of the boiler to model all the dynamics of the boiler. These states are Drum Pressure, Total Volume of Water, Steam Quality, and Steam Volume in Drum. The state of drum pressure is also an output. The other output is the drum level, which is modeled as an independent nonlinear function of all the states. The overall model [33] was derived using laws of conservation of energy, mass, and momentum for individual components as well as for the whole system globally. The governing equations of the system take the following nonlinear state-space form: where the coefficients a ij s are where P, V wt , α st and V sd are the state variables. Other variables that are present in the model are computed using the following relations: In model (1), the output level is not explicitly defined as a state. One way would be to use Lie derivatives approach and represent level as one of the states in a new model, but it is not doable given the complexity of nonlinear equations. To make it easier, the level is defined as a separate output variable which is computed using the following set of equations: where LUT i 's are lookup tables, which for boilers, are referred to as steam tables that provide values of ρ i 's and h i 's as a function of the sate of steam pressure P. A polynomial function of a second or third-order is sufficient to interpolate the values of the steam table.
The given model (1) assumes a typical state-space form, but the real complexity lies in the algebraic manipulation of nonlinear coefficients of (2). One issue in the MIMO system like the boiler is that there exists a mathematical coupling between all the inputs and outputs of the system. So it is necessary to determine which input-output pair has to be formed for constructing the control loops. This can be done by analyzing the mathematical model or doing hit-and-trial simulation to analyze the impact of each input on all the outputs. For the case of boiler, it is observed that the following input-output pairs are most efficient in executing a controlled response of boilers: • Feedwater rate -level control. • Heating rate -pressure control. The level control is implemented under the name of threeterm control, in which the feedforward element of the steam flow rate is included to nullify the variations in the steam variable. Figure 1 shows the complete control loop structure. The feedwater rate is used in a feedback control loop to stabilize and regulate the level output against dynamic disturbances of steam flow rate: where ''f L '' is the control function to be formulated in the next section. Similarly, the pressure control loop is designed using the control input of the heating rate: where ''f Q '' is the control function to be formulated in the next section. Given the complex nature of the model, as evident from (1)- (15), analytical derivation of any kind of control technique is impossible. Typical analysis of Lipschitz properties of the functions is also not doable, and deriving an analytical solution of the system response is barely possible due to such nonlinear structure of the model.

IV. NUMERICAL FORMULATION OF SMC FOR BOILER
The basic principle of SMC design is to transform an 'n dimensional tracking problem into a one-dimensional stabilization problem. This achieves a controlled response in two phases, commonly known as the reaching phase and the sliding phase. Various methods are available that proposes different designing schemes for both the manifold design and the control law. For this work, we use Utkin's approach for the basic derivation of 1-SMC. However, it is to be noted that most of SMC design methods are based on the models having strict feedback integrating form (normal form or canonical form). For such models, the relative degree of output variables is at least two, hence designing a stable sliding manifold is simple. In the case of the boiler model, all the outputs and states have a relative degree of one, and hence we use a trivial formulation of sliding manifold using one variable. Later, we will upgrade the sliding manifold to include the two output variables. At this point, we avoid designing manifold as a linear combination of all states, because for highly complicated nonlinear models, it is unclear whether the linear combination of all the states would produce a stable sliding manifold or not. That's why we restrict the design process for the two variables, which are in fact the outputs of our system.
Consider model (1) of our system in the following form: whereṁ df =ṁ fw −ṁ s is differential feedwater rate, and a ij 's are functions of all states as given in (2). V 0 sd is taken 0 for simplicity, however in real operation it has some value.
The above model takes the following notational form: where, VOLUME 8, 2020 For ease of notation, we omit the index 'x', and use the following form: The state (19) can be rewritten as: The equation (20) uses the inverse of matrix E, which is only calculable in the numerical form given the nonlinear structure of matrix E. Such modeling form of boiler (20) is of great essence in the formulation of numerical SMC. Since it is impossible to have E −1 and (20) in analytical form, the only way to implement this equation is in the numerical way, i.e., by computing (20) numerically at each instant of a simulation or sampling instant of a real-time operation. So, at each instant, we produce the numerical form of the equation (20) and do the subsequent computations of numerical SMC design. There are two outputs to be controlled, the drum level 'L , and the steam pressure 'P . For now, we focus on control of the pressure. The level would be controlled by using a PID with nominally tuned parameters.
The pressure is the second state in the state vector, and so is given as: where, C p = [0100]. For our system, the input relative degree of all the states and outputs is one. In this respect, we proceed with a trivial form of the sliding manifold, where the sliding variable 's' is simply a function of the error between output and setpoint. This form of the sliding manifold is also used in [6], [37], [38], and can be considered as a sliding point instead of a manifold and is commonly used in Lyapunov control design. We begin with this sliding point as a bottom-up approach for more complex manifold design later in this article. Consider the sliding variable as a function of output pressure.
whereP = P − P d , and P d = 8500 is the desired state value of steam pressure. We break the input vector into two inputs as follows: In this way, the 'B' matrix is dissociated as: For now, we consider the caseṖ d = 0, i.e. constant desired pressure and so, the sliding variable derivative is given as: Computing u eq using Utkin's approach [1], Decomposing input 'u' into u Q and u SF .
Since the modeling inaccuracies cannot realize an ideal equivalent control of (29), the dynamics of the system get partially canceled out and remain there in the form of a perturbation. This perturbation factor is taken into account while designing the switching gain 'k .
For equivalent control, we use the nominal form of the model (19): where, E o , A o and B o are nominal matrices that are dependent on the model parameters selected, for example, by using average values of the uncertainty bounds. Since the actual parameters of the system cannot be known, the differences between the nominal matrices and actual matrices will always be present, but the good point is that such differences are bounded and numerically calculable. We compute these bounds in Section VI. In order to use Utkin's method of SMC [1], we need to determine the switching gain 'k of (29). We can do so by using the nominal form of (29) and (31) in the derivative of sliding variableṡ of (27). Let, In order for the sliding surface to remain attractive, a quadratic Lyapunov function of the following form can be used to ensure the reachability of the manifold and hence the stability of the overall system.
which is only possible if the following inequality is satisfied: In this work, the gain 'k' is determined numerically at each instant, using the following equation.
where k 0 = |γ | and η is arbitrary that determines the reaching time of error trajectories and can be manipulated at the expense of control energy. Equation (39) constitutes a numerical form of SMC for pressure output, where all the terms of expressions are available in numerical form and so is k o . The parameter k o is in fact a function of the states i.e. k o = k o (x) , and change dynamically with all the states. Hence the parameter k = k 0 + η can be regarded as an adaptive gain of the switching control of the SMC.

A. SMC FOR LEVEL CONTROL
The numerical design of SMC for the level output uses a modified procedure of the SMC design for the pressure output. This is because the output level is not explicitly defined as a state as the output pressure is. The level is modeled using the nonlinear equation (8), which further depends on expressions (2) and (9). Combining all the expressions analytically to produce a single equation is not possible, so we use a general expression to denote level output as: where 'x' is a state vector, and 'h' is a composite function that includes all the functional dependence of level on all the states.
Here, we refer to the same sequence of SMC derivation of the previous section where equations (28), (29) and (38) are the key components of the numerical formulation of SMC. The derivation of SMC, according to these equations, is contingent on the format of the pressure output, which is in linear form i.e. y = P = C p x. Hence it is required to convert (40) into a similar linearized form, so that an equivalent of equations (28), (29) and (38) can be derived for output level as well. For that, we first examine (25) and write it in a general form: whereĀ,B,C can be any general matrices, but for our casē The variable 'y is the concerned output, which is either level or pressure or both as a vector. From (40), the time derivative of output level takes the following form.L = ∂h L ∂xẋ where ∂h L /∂x is not available in an analytical form due to the complexity of the nonlinear function h L . Here we use the approach of numerical differentiation using forward difference method, where ∂h L /∂x is numerically represented in the following form.
Essentially (42) represents the gradient of h L using the numerical approach. The gradient is computed in the form of an algorithm rather than a mathematical formulation. We use the Ab Initio method (or known as the First Principle) of differentiation to form this algorithm, which can be generalized into vector form or even a matrix form for computing vector and matrix time derivatives. The algorithm is, now onwards, termed as Ab Initio algorithm and is as follows: • Slightly increment each state, i.e. x i = x i + x • Measure the corresponding change in L, i.e. L. • Compute L/ x i . • Repeat the above steps for all states and create the gradient vector: =Ĉ ∇L • The numerical value of ∇L, when multiplied withẋ gives an estimate ofL. VOLUME 8, 2020 Using the above algorithm gives us a linear formL =Ĉ ∇Lẋ similar toṖ = C pẋ , which can be analytically used to represent the level in the form of (41): Since we useṁ df to control level, we dissociate input vector in the following form to have an independent representation of differential feedwater rate U dF =ṁ df .
The equation (44) is applied onwards in the SMC design for the level output. A simple test was conducted to prove the applicability of this method. A disturbance signal of the sinusoid is applied using the input of steam flow rate at time t = 200, and derivatives of level are recorded using the two methods; using the Ab Initio algorithm with (43) and using the backward-difference method to computeL directly. Figure 2 shows that both methods produce almost the same results proving that computingL using Ab-Initio algorithm works accurately. We use this derivative of the level to produce the equivalent of equations (28), (29) and (38) which takes the following form: The above numerical formulation of SMC based level control forms the basis of generalized numerical derivation of SMCs for the systems whose output is a nonlinear function of all the states. The fundamental idea is that we need to derive a linearized form of the output similar to (43) so that it can be used easily with the available methods of SMC algorithms. We generalize this method in detail in Section VII.

V. SECOND ORDER SMC
Higher-order sliding modes are usually designed with the aim of improving the SMC control on the following grounds: • Chattering reduction.
• Reduction of the sliding mode order.
• Finite-Time convergence. Initially, HOSMCs were introduced as Second-Order SMCs for which many authentic algorithms have been introduced so far, each having its own significance. Essentially, in 2-SMCs, the error trajectories are required to satisfy the following sliding constraint: The relative degree of sliding variable s can be both 1 or 2. In the case of a relative degree of one, the constraint is satisfied using the derivative of the input instead of regular input. Since our case system (boiler) involves output variables of relative degree one, we focus on the 2-SMC formulation for the sliding variable of degree 1. Consider the following expression of the sliding variable: Taking first and second-order derivatives of s, which can be rewritten in notational form as: where the following bounds hold: For the case of the boiler, we incorporate the system equations into the above equations to determine the time derivatives of 's . Consider a sliding variable 's defined as a linear combination of outputs: For now, we derive 2-SMC for output pressure y = P. The level will be controlled using a secondary controller (for example, PID).
Taking 1 st and 2 nd order time derivatives and incorporating (20): whereu SF should be zero, as the concerned input is u Q , while other inputs are considered constants. The above equation is in the following form: where, Equations (59) and (60) are sufficient for the formulation of 2-SMC. Particularly, these equations can be used to numerically compute the bounds of variables α(t, x, u) and β (t, x, u). Once the bounds are available, any form of 2-SMC can be designed and implemented. In this work, we use the STC algorithm to prove the compatibility of our numerical approach with the 2-SMC. The STC law is formulated as: where the following constraints should be satisfied:

VI. NUMERICAL COMPUTATION OF THE BOUNDS AND THE TIME DERIVATIVES
This section discusses a very important aspect of this work, which is the numerical computation of bounds and time derivatives that are involved in all SMC formulations.
From the above equations, especially (36), (59) and (60), it is evident that we may not need the actual values of γ (x, t) , α (x, t) and β(x, t), however the bounds of these variables are important in the SMC design process. Once the bounds are available, the constraints (38) can be easily fulfilled to implement first-order SMC or HOSMC algorithms. The bounds are not known apriori and cannot be computed analytically for complicated models like that of boiler. Hence, we use a numerical approach to compute these bounds. For our method, we rely on the assumption that the parametric uncertainties in the system's original parameters are all known. These uncertainties are also bounded with known constraints in the following form: where 'κ is the vector of systems static parameters that totally depend on the particular structure and configuration of the system. For boilers, each element κ i of κ belongs to the set , i.e.
The parameter's nomenclature can be referred to in the Appendix IX-A. The parametric uncertainty (63) can be written for all the individual parameters as: These uncertainties are transformed into the bounds (38) and (53). The transformation is not straightforward and cannot be done analytically, and hence a numerical computation of bounds is carried out. For this, we evaluate γ (x, t) , α (x, t) and β(x, t) at a given operating point using a range of system models chosen from a discretized set satisfying the given domain (65). We keep track of how the maximum and minimum values of γ (x, t) ,α (x, t) and β(x, t) vary with respect to all extreme values of the uncertainties. The maximum of these variations is termed as deviation factors ( γ , α, β) that are also used to make a robust estimate of the bounds. This approach can be stated using the following algorithm: 1) Simulate the system in 'N' number of iterations, where 'N' is the maximum number of combinations of κ l and κ m of (65). Mathematically N = 2 N p , N p being total number of uncertain parameters. We may discretize the set (65) into more points at the expense of computational effort to produce more accurate bounds. 2) In each simulation, iteratively set κ i to one combination of maximum and minimum bounds (κ l and κ m ) of (65).
b. For 2-SMC, use: During regular operation, we adaptively compute the bounds k for 1-SMC, and m , M , for 2-SMC. For example, at each time instant, we compute the numerical value of β (t, x, u). We adaptively set m as the lowest value of β (t, x, u), and M , as the highest value of β (t, x, u) plus/minus the deviation factors β l or β m that arises due to parametric uncertainties of the model. An example Matlab implementation would be in the following way: If In a similar fashion, we numerically compute the bounds k and of γ (x, t) and β(x, t), respectively.

A. COMPUTING TIME DERIVATIVES OF MATRICES
The expressions, (59) and (60), involve time derivatives of matrices in the formĖ −1 ,Ȧ,Ḃ SF ,Ḃ Q . We compute these derivates numerically, at each instant, using the same Ab Initio algorithm that we used forL. Since,L is a scalar, the approach needs to be applied for all the individual elements of the matrices, or we can do a matrix implementation as well with a little increase in the complexity of the programming code. However, using a matrix approach results in a three-dimensional cube comprising four matrices due to the four states of the boiler system. The cube is multiplied withẋ, which gives time derivative of the matricesĖ −1 ,Ȧ,Ḃ SF ,Ḃ Q .

VII. GENERAL FORMULATION OF NUMERICAL SMC
After developing SMC for the boiler nonlinear model, we can generalize the approach for nonlinear models of any form.
Consider a general form of a nonlinear model: where x ∈ R n , u ∈ R m ,u d ∈ R q and y ∈ R p . The input u d is assumed to be a measurable disturbance.
Let the corresponding nominal model in the following form:ẋ Consider, a sliding manifold defined as a linear combination of the error states: where we have assumed that the sliding manifold Gx is inherently stable. Such stability must be ensured based on the knowledge of the system or by-hit-and-trial simulations.
The relative degree of 's is assumed to be 1. The SMC formulation for the case of the relative degree of 2 is easier and is not discussed in this work. Also, we only consider constant setpoints, i.e.ẋ d = 0. However, the approach can be extended to the cases with variable setpoints with little modification. Differentiating the sliding variable w.r.t time, The above sliding manifold can be easily transformed as a linear combination of output errors by using the following relation:ṡ where G = GĈ ∇y . The vectorĈ ∇y is computed using the same Ab-initio algorithm that was used to compute the gradient of level ∇L ij or matrix element ∇e ij in the previous section. Using (68) in (71) Using, u = u eq − ksign (s). Also, dropping the index (x,t), we can rewriteṡ equation as: where the input u, has been incorporated from the nominal state (69). Let, To maintain reachability condition of the sliding manifold, the following condition must be fulfilled: A trivial fulfillment of the above constraint can be carried out using the following: where and η is arbitrary.
Similarly, for second-order SMC, we need to compute numerical bounds of variables (α, β) in the equation ofs. Suppose that the sliding manifold is a Hurwitz stable linear combination of outputs i.e., s = Gỹ, then using (68), we can computes in numerical form.
Taking the second-time derivative: where,ẋ is available from (68). The matrices f x , g x , and φ x are the Jacobians of the functions f , g and φ respectively with respect to all the states. The Ab-initio algorithm of Section VI-A is used for their computation, which because of brevity, is not discussed again. The numerical computation of G x is also carried out using a similar numerical approach that is presented in detail in Section IX-B of the Appendix.
Rearranging (85), Comparing above with the general forms =ᾱ (t, x, u) + β (t, x, u)u(t), we have: The algorithm of Section VI can be generalized as follows: a. For 1-SMC, usek where γ is numerically available from (77). γ m is the deviation factor numerically computable using the algorithm of Section VI. b. For 2-SMC, use: whereᾱ andβ are numerically available using (87) and (88). β m , β l and α m are the deviation factors numerically computable using the algorithm of Section VI.

VIII. SIMULATION RESULTS
We simulate all the above-presented controllers while investigating the responses of the two outputs of the boiler. The input of steam 'ṁ s ' is used as the main disturbance agent in the boiler system as any variation in the steam can cause all the unstable boiler dynamics to go unbound. We change the steam flow in the form of a ramp at t = 100s. All control implementations are individually discussed in subsequent subsections. In order to create uncertainties, we randomly change the parameters of the boiler by ±10%. For simulation, we run two models of boiler simultaneously; one represents the actual system whose parameters are changing by ±10%, the other represents the nominal model with constant parameters whose equations are used to determine SMC equivalent control law. The aim is to stabilize pressure at 8500 and level at 1. The values of all boiler parameters are adopted from [33].
A brief overview of how the controllers are set is presented hereafter.

A. 1-SMC FOR PRESSURE CONTROL
We have first used 1-SMC to control one output to observe how the proposed approach is working. In this regard, we take a sliding manifold of (22) to control pressure using the input of the heating rate 'Q . The level is controlled by a nominal control (PID) via an input of differential feedwater rate 'ṁ dF '. Figure 3 and Figure 4 shows the simulation results with   The 1-SMC approach can be upgraded to sliding manifolds of various forms. In order to tackle out the disturbances, Integral SMC (ISMC) is a widely recognized SMC formulation, in which the order of sliding variable is increased by incorporating an integral function of error. In this respect, the ISMC manifold takes the following form: Using this form of sliding variable 's , the mathematics of the numerical design of SMC is not affected and hence the same program is re-simulated with the upgraded sliding manifold. Figure 5 shows the simulation response under the same operating conditions, where the following parameters were considered in order to design the manifold: The simulation again proves that the numerical approach of designing integral SMC is working perfectly while stabilizing all the concerned variables at respective setpoints without violating safety constraints.

C. 1-SMC FOR LEVEL CONTROL
In the next phase, we switch the control loops. We use SMC to control the output level while using a nominal PID for pressure control. For level control, we test another formulation of SMC, which is described in Khalil's text [39].
where, µ (x) = µ 0 + ρ (x) + ρ Where ρ is the same deviation factor that is determined as ρ = 5.8059 using the algorithm of Section VI. µ o is arbitrary and has the same role as η, i.e., to determine the convergence time of error trajectories towards the sliding manifold. The parameter ρ(x) is derived from (44) such that the following is satisfied.
This control formulation also works accurately. Figure 6 shows the simulation result where the output level was controlled at the desired setpoint of L d = 1.

D. COMBINED 1-SMC FOR PRESSURE AND LEVEL
We implement MIMO control using sliding manifold as a linear combination of output errors. The sliding gains in multi-variable manifold should be selected such that the resultant sliding dynamics remain stable in the sliding phase. This is a challenging task for a complex nonlinear model, where the model is not available in normal or canonical form. Hence, the only way is to randomly select the gains of the manifold and simulate the results to see the stability of manifold for the selected gains.
In MIMO SMC formulation, the output pressure is controlled by heating rate Q in one SMC, and drum level is controlled by differential feedwater rateṁ dF using another SMC. To realize this, a 2-dimensional sliding surface is formulated. Here we point out that only 2 outputs have been used for the boiler in the manifold design, where a more generalized approach would be to use all the states. This is not possible for the boiler case because, due to the nonlinear nature of the model, the x id s (the desired setpoints of all the individual states corresponding to output setpoints) need to be determined for every simulation instant, which will be analytically impossible and even computationally complex using the numerical approach. First, we investigate how the two SMCs are independently designed without any coupling between the controls. This can be done by enforcing the de-coupling condition as g P2 = g L1 = 0. Figure 7 shows the simulation results.
In the next phase, we relax the decoupling condition by choosing g ij s = 1, ∀i, j which turns out to be a stable manifold. Figure 8 shows the simulation result of this approach. The offset in level output can be trivially neutralized using an integral state x L i = t 0L (s) ds.

E. ANALYSIS OFṠ ANDS NUMERICAL COMPUTATION
As a validation step, we test a simple model of second-order in order to analyze how the numerical computations ofṡ ands compare with their respective analytical computations. We first derive the expressions ofṡ and s analytically, and then validate the accuracy of the numerical approach for the equations (85) and (72). The model selected is from the text [40], which is as follows.ẋ The following SMC control of the example 1.3.2 from [40] has been used to track r = sin(t).
Let's define the output variable as a nonlinear function of states.
We use a sliding variable as a function of output: Taking first and second derivatives while considering y d to be constant.

F. SECOND ORDER SLIDING MODE CONTROL
We simulate Super Twisting SMC as presented in Section V using the proposed numerical approach. Only the output of pressure is controlled, while the output of the level is stabilized using a nominal PID. The mathematical form of Super Twisting SMC control law, as given in (61) is used to simulate the results. To satisfy the constraint equations (62), we use the following equations: To satisfy W > m , we select W = m + 1. To satisfy λ 2 ≥  Figure 10 shows the simulation results. It is evident that the Super Twisting SMC response VOLUME 8, 2020 shows a diminished chattering as compared to all other forms of SMCs.

G. BRIEF DISCUSSION
The simulations (Figure 3-Figure 10) show that the introduction of change in steam dynamically affects all the variables; however, they remain within safety limits without violating any physical constraint. Also, it is evident that the SMC successfully stabilizes both outputs, the pressure, and the level, which after undergoing temporary dynamics, maintain their respective setpoints. All the simulations show that the inputs are able to regulate the concerned outputs at their respective setpoints. The inputs of heating rate and feedwater rate change dynamically for a small interval until they attain new equilibrium points, which neutralizes the effect of steam disturbance while minimally affecting the output variables. Both inputs follow the ramp like the trend of steam flow, which is intuitive from the model (1), as the dynamics of these inputs have to change equally to keep the right-hand side of (1) close to zero. With an increase of 45% in steam flow, the variations in level output are confined within a bound of 5% while variations in pressure are confined within a bound of ±1 kPa. (<0.01%) in all the simulations.
As an extension of this work, we may integrate an observer with the SMC and implement the proposed approach to carry out a combined observer and controller design. It may also be considered to apply the proposed approach with other non-linear control techniques since the numerical design approach is not limited to SMCs. Particularly, back-stepping control design would be a good candidate to apply the numerical approach for system models of complex form.

IX. CONCLUSION
Analytical formulation of Sliding Mode Control (SMC) with complicated dynamic equations is not possible. This is because algebraic manipulation of state equations becomes undoable if the nonlinearities complexity exceeds a certain level. Given a system's model in affine form, a numerical approach can be used to design the SMC control law. The control needs to be computed at each instant of simulation by applying a numerical algorithm. The proposed numerical algorithm of this work solves the SMC design problem by doing two important roles; computing time derivates of nonlinear functions using approximate gradients, and computing bounds of variables that are involved in time derivatives of the sliding variable. For the computation of the approximate gradients, a numerical approach (Ab-initio algorithm) based on forward-difference is introduced. For the computation of bounds, the parametric uncertainties bound of the system being known apriori, are numerically transformed into bounds of the variables of sliding variable derivatives. Also, a deviation factor is incorporated into the bounds that take into account the parametric uncertainties of the system. Once the time derivatives are computed, and all the bounds are numerically available, it remains no more difficult to derive and implement any form of SMC controllers. The approach has been corroborated with the application of few forms of First-order SMC, as well as Super-Twisting Second-order SMC. The simulation results, with a case study of the boiler, prove the applicability of the proposed numerical approach.

A. NOMENCLATURE OF SYMBOLS AND VARIABLES
See Table 1.

B. DOUBLE DERIVATIVE OF FUNCTION h(x, t )
Consider a case in which the aim is to numerically compute the gradient of a vectorĈ ∇y .
where G is a constant matrix.Ĉ ∇y itself is regarded as the gradient of h along the state vector x.
The further gradient of the above vectorĈ ∇y , i.e.
Which is evaluated using the Ab-initio algorithm by using minute increments to numerically compute the individual elements of the above matrix. The double derivatives of elements are computed as follows.
Consider an (i, j) element of the above matrix.
The point (a, b) represents the current value of (x i , x j ) at a certain instant. To compute the double derivative of 'f', we apply very small increments, a → 0, b → 0, in (x i , x j ) while keeping the rest of the states constant, shown at the bottom of the page. In the case of (i =j), i.e., for the diagonal elements, the above gets simplified.
where h → 0, and is assigned a very small numerical value in the simulation. The rest of the states are kept constants.
Using the above two equations, all the elements of ∂(Ĉ ∇y ) ∂x are numerically computed and used thereafter.