Model Predictive Control for First-Order Hyperbolic System Based on Quasi-Shannon Wavelet Basis

: In this paper, we consider a class of ﬁrst-order hyperbolic distributed parameter systems. Our focus is on the design of a new class of model predictive control schemes using a quasi-Shannon wavelet basis. First, the ﬁrst-order hyperbolic distributed parameter system is transformed into an equivalent system using collocation techniques for the approximation of spatial derivatives and Euler forward difference method for the approximation of the time component. Then, a model reduction method is applied to obtain a reduced-order system on which a nonlinear model predictive controller is designed through solving a nonlinear quadratic programming problem with input constraints. For illustration, the temperature control of a ﬂow-control long-duct heating system is considered to be an example. A comparative simulation study is conducted to demonstrate the effectiveness of the proposed method.


Introduction
Modeling and control of systems governed by partial differential equations (PDEs) have attracted the interest of many researchers. See, for example, [1] and the references cited therein. These systems are known as distributed parameter systems (DPSs). Among the controller synthesis methods for DPSs, there are direct design strategies [2], post-centralized strategies [3], and prior-centralized strategies. For these three strategies, the prior-centralized strategy is easier to be understood, so it has been widely used to deal with DPSs. With the prior-centralized approach, a DPS is first approximated by a system of ordinary differential equations, referred to as a lumped parameter system (LPS) [4]. However, the obtained LPS tends to have a high dimension, and hence it is essential to get a low-order representation of the LPS and use it for the design of the controller. There are many effective control schemes in the literature for the design of the controller.
For a system described by hyperbolic PDEs, the eigenspectrums of its differential operator contain almost all the energy. As a consequence, the order of the model is higher than that of the parabolic system [4]. The classical finite difference method (FDM) has long been recognized as a popular numerical method for DPSs among the control community. For example, several FDM numerical models have been developed to deal with various control problems, including optimal control problems governed by first-order hyperbolic systems. A coupled system involving PDE and ODE is transformed into a zero-order neutral system with distributed delays [5]. An optimal boundary control problem, where its underlying dynamical system is governed by linear hyperbolic partial differential equations with a turnpike structure, is being applied to solve a control problem in gas pipeline operation [6]. A boundary input-to-state stability of a class of coupled systems involving hyperbolic PDE and ODE in the presence of uncertainties and external disturbances is investigated and the results obtained is applied to an interconnected continuous stirred tank reactor coupled with a plug-flow reactor [7]. Furthermore, the weighted residual methods (WRMs), such as the Galerkin method, collocation method, and finite element method (FEM), have also attracted the attention of many researchers among the control community. WRM has been widely used in various industrial applications due to its high computational efficiency and flexibility in handling complex geometries and boundaries [8]. For control problems governed by various DPSs, the weighted residual method based on a wavelet basis is found to be effective, showing great potential [9]. The collocation method involves finding a functional basis for the solution space of the equation, then making the residual zero at n points within the domain [10]. Generally, wavelet-collocation method (WCM), created by choosing a wavelet and some grid structure, is most appropriate for solving nonlinear PDEs with general boundary conditions [11]. The computational model of hyperbolic wave equations is proposed by using Haar wavelet operation matrix [12]. A non-Fourier heat conduction problem is solved through the use of Haar wavelet basis implicitly in combination of finite difference and pseudospectral methods, thereby avoiding numerical oscillations near wavefront jumping discontinuities [13]. Although the Haar wavelet function has an explicit expression, it is non-differentiable. Thus, an integral smoothing calculation is introduced. Furthermore, since the Haar wavelet has only one vanishing moment, it does not perform well for the approximation of a smooth function. The Pennes bioheat transfer equation is solved by the WCM and a multiresolution analysis is conducted using Daubechies wavelets [14]. The Daubechies wavelet is tightly supported, smooth, and has n-order vanishing moments, but its differential operator is difficult to obtain due to the lack of analytical expression. By combining the operation matrix of integration for Legendre and Chebyshev wavelets with the collocation method, an approximate solution of a one-dimensional hyperbolic PDEs with given initial conditions is derived [15]. By the operation matrix of fractional integration combined with the Galerkin method and the collocation method, a coupled system of Burgers' equations with the time-fractional derivative of Gegenbauer wavelets is obtained [16]. The operation matrix of Hermite wavelet is derived, and on this basis a direct numerical method is developed for solving optimal control problems [17]. By combining WRM with Shannon wavelet basis, a numerical method is proposed to calculate the stationary solution of Fokker-Planck-Kolmogorov equation [18]. Shannon wavelet is continuous, smooth, and infinitely differentiable, but it is not tightly supported. For the quasi-Shannon wavelet [19], it is also tightly supported.
Model predictive control (MPC) is a successful method for the control of industrial processes modeled as LPSs or DPSs with constraints and time delays. For the control of hyperbolic systems, such as long duct heating systems [20], and convection-dominated chemical reactor [21], the need to handle constraints and time delays is unavoidable. The key step in the design of an MPC algorithm for DPSs is to build a prediction model through model reduction. For DPS, a linear MPC scheme using Cayley-Tustin transformation is proposed in [22]. The FDM has also become a popular method for constructing prediction models. MPC using FDM's prediction model has been successfully applied to several industrial processes, such as convection-dominated transport-reaction processes [21], cocurrent fluid flow problems [22], non-isothermal plug-flow reactor [23,24], wind tunnel propagation delay problems [25], and open irrigation canals [26]. The characteristic line method (CLM) is a spatio-temporal discretization method. It has been applied to the study of control problems, where the SO 2 concentration is to be controlled in a circulating fluidized bed [27]. It has also been applied to uncertain hyperbolic systems [28], and a counterflow plug-flow reactor [29]. On the other hand, model predictive control based on the weighted residual method is being developed for hyperbolic systems.
The control problem of selective catalytic reduction in diesel-powered vehicles is solved by CLM combined with spectral decomposition [30]. Wavelet-based numerical methods have great potential in MPC application design, but there is little related literature. A prediction model based on Haar wavelet basis has been applied to the design of the MPC controller of discrete DPS [31]. However, in order to improve the computational efficiency, it is necessary to develop a more suitable wavelet-based method.
In this paper, the WCM is used to construct the reduced-order model of first-order hyperbolic system. The interval quasi-Shannon wavelet scaling function is used as a basis function to suppress the boundary error. Then, a nonlinear MPC scheme is developed based on the above prediction model. Nonlinear MPC schemes are also developed based on the models obtained by FDM and the Haar wavelet transformation proposed in [31]. They are applied to solve the disturbance rejection problem of long duct heating systems. Their computational efficiency has been compared.
The rest of the paper is organized as follows. In Section 2, the first-order hyperbolic system and its control objective are introduced. Then, the interval quasi-Shannon wavelet is introduced, and a nonlinear model predictive control scheme is developed based on the wavelet-collocation method for the first-order hyperbolic system in Section 3. In Section 4, a simulation study involving a long duct heating system is performed to verify the effectiveness of the algorithm proposed in Section 3. Some concluding remarks are given in Section 5.

Problem Formulation
We consider the following class of distributed parameter systems described by first-order hyperbolic evolution equations: subject to initial condition y(x, 0) = y 0 (x) and boundary condition y(α, t) = L(t), where x ∈ [α, β] ⊂ R is the spatial variable, R is the set of real numbers, t ∈ [0, ∞) is the time variable, y(x, t) ∈ H [(α, β), R] is the controlled variable, H [(α, β), R] is an infinite-dimensional Hilbert space containing integrable spatial differentiable functions defined in the interval [α, β], y(x, t) is differentiable with respect to x and t, y t (x, t) = ∂y(x, t)/∂t, y x (x, t) = ∂y(x, t)/∂x, u(x, t) is a smooth and bounded manipulatable variable, f 1 , f 2 are sufficiently smooth nonlinear functions of x and t, L(t) is a smooth and bounded function of time t, f 1 2 ≤ A < +∞, f 2 2 ≤ B < +∞, where · 2 is L 2 -norm. The optimal control problem may now be stated as follows: where J(u) is the objective function, and r(x, t),ŷ(x, t), u(x, t) and ∆u(x, t) are, respectively, the reference trajectory, the prediction output and the control increment at position x and time t. u min and u max represent, respectively, the upper bound and the lower bound of the control variables, while ∆u min and ∆u max represent, respectively, the upper bound and the lower bound of the control increment. λ is the weighting coefficient.
To proceed further, let us consider an actual situation, which is referred to as Problem (P), as stated below.
Consider the situation where there are M = 2 J + 1 sensors, which are uniformly located, and l actuators, and let u( , where H is the standard Heaviside function.
The discrete forms of the system and the cost function will be constructed using the wavelet-collocation method to be described in Section 3.

Model Predictive Control Using Quasi-Shannon Interpolating Wavelet
This section describes the principle and design process of the proposed MPC schemes using a quasi-Shannon wavelet basis. After a brief introduction of the principle of FDM and WCM, the given first-order hyperbolic system is transformed into an equivalent system using a collocation technique for spatial derivatives combined with the Euler forward difference method for the approximation of the time component. On this basis, the resulting prediction model is constructed. Since the quasi-Shannon wavelet scaling function has advantages in approximating continuous functions, the quasi-Shannon wavelet scaling function on a finite interval is selected as the basis function. On this basis, a prediction model is constructed, which is then applied to design the MPC scheme for the first-order hyperbolic system such that the discrete objective function is minimized.

Model Reduction for First-Order Hyperbolic Systems
There are various model reduction methods in the literature, such as the eigenfunction method, Green's function method, and WRMs [4]. For WRMs, they include Galerkin method, the collocation method, the FEM, and the spectral method. We will consider the latter two methods because the eigenfunction method and Green's function method are not suitable for nonlinear DPSs.

Finite Difference Method
The FDM is the earliest and most popular approach to approximate PDE by difference equations (in which the derivatives are approximated by finite differences) [32]. Based on the Taylor series expansion, the PDE system is approximated by a system of difference equations. These difference equations are obtained based on forward, backward or central differences. The accuracy of the algorithm is determined by the grid size and the scheme for constructing the difference equations. This scheme must satisfy the Courant-Friedrichs-Lewy (CFL) condition to ensure the stability of the discretization.
We consider the first-order hyperbolic system (1). The mesh of the solution region is as shown in Figure 1. Let the space discretization interval be ∆x and the time discretization interval be ∆t. Then, x m = α + m · ∆x, m = 0, 2, · · · , M − 1, t n = n · ∆t, n = 0, 1, · · · , ∞, and y n m = y(x m , t n ), f | n m = f (x m , t n ) are, respectively, the approximations of y(x, t) and f (x, t) at x m and t n . The upwind discretized form of system (1) can be expressed as: where the accuracy of the upwind discretized form is first-order, and the CFL condition is ∆t ≤ ∆x/max{ f 1 | n m }. The FDM-based MPC scheme will be compared with the proposed method in Section 4.
FDM has two disadvantages. It is more complicated when solving PDEs on irregular domains, and the mathematical analysis is difficult for nonlinear PDEs with FDM [11].

Spatio-temporal reconstruction
Step length Meshing Difference format Spatio-temporal discretization using difference Equations solving by elimination or iteration Method Distributed parameter systems

Spatio-temporal reconstruction
Step length Meshing Difference format Spatio-temporal discretization using difference Equations solving by elimination or iteration Method

Wavelet-Collocation Method
The wavelet-collocation method is derived based on the WRM, which is one of the most effective numerical methods for solving PDE systems. The numerical method obtained based on WRM is conceptually different from those obtained using FDM because WRM assumes that the solution can be expressed analytically. It is well known that a continuous spatio-temporal variable can be approximated in terms of a set of spatial basis functions Equation (4) is truncated as: For (1), we can express the residual R(x, t) as: To minimize R(x, t), we let where {ϕ i } M−1 i=0 is a set of weighting functions. The residual R(x, t) is converted to its projection onto the space of weighting functions. The key to the accuracy and efficiency of WRM is the selection of the basis functions and the weighting functions.
The two most popular approaches to procure the weighting functions are Galerkin method and the collocation method. can be classified as global and local types. Typical basis functions include Fourier functions, orthogonal polynomials, and wavelet functions. Various model reduction approaches have been proposed through the application of Galerkin method and the collocation method using different basis functions [4].
Because of the collocation nature of the WCM, it is a simple task to deal with nonlinearities. Compare with wavelet-Galerkin method, the collocation method is more practical for implementation [11] (see Figure 2).

Interval Quasi-Shannon Wavelet
In this section, the quasi-Shannon wavelets are used as the basis functions of the proposed prediction model taking into account its normalization conditions. Then, the choice of the parameters is discussed. To eliminate the boundary effect, the interval quasi-Shannon wavelets are adopted.
Shannon wavelet functions are, which have analytic expressions, orthogonal, smooth, and differentiable. However, they are not tightly supported. To derive a wavelet function with compact support, Shannon wavelet is multiplied by Gauss window function, forming a semi-orthogonal wavelet function with compact support. This wavelet function is called the quasi-Shannon wavelet [19]. The scaling function φ ∆,σ and the wavelet function ψ ∆,σ of the quasi-Shannon wavelet are, respectively, given by where ∆ is the spatial mesh size, while σ is the window size. The first-order derivative of the quasi-Shannon scaling function is depicted in Figure 3 obtained by Matlab 2016b, and it is expressed as: To guarantee normalization conditions being satisfied, the following lemma [19] is given.
The choice of the parameterr and the nature of the curve being described with wavelet as the basis function will affect the accuracy of wavelet representation. A guide for the choice ofr and the resulting mathematical estimation of the approximation error are given in [19]. Whenr → ∞, the tight support property of the wavelet is completely lost. On the other hand, whenr → √ 2/π, the normalization conditions will fail to be satisfied.
Empirically, the parameterr should be chosen in the interval [2.1 2.3] (respectively, [2.6 3.8]) when approaching a linear function (respectively, a nonlinear function). The wavelet transform is defined in the bilateral infinite interval, but the signals appearing in real life are usually in a finite interval. The interval wavelets are often used to reduce boundary effects. The interval quasi-Shannon wavelets based on the difference method [33] are given by where j is an integer, φ(2 j x − k) is the quasi-Shannon wavelet scaling function, N is the number of the outer wavelet scale function on each side of the boundary.

The Prediction Model Developed Based on WCM
Based on WCM under the situation referred to as Problem (P). The spatio-temporal variable y(x, t) in (1) is approximated by where y M (x, t) is the approximation at the spatial location x and the time t. Let φ(x) be the scaling function of the quasi-Shannon wavelet, and let φ j | x−x m = φ j (x − x m ). Then, the first derivative of the spatio-temporal variable y(x, t) at any time t is given by By Euler forward difference method, it gives explicit format and a first-order accuracy approximation to y t (x m , t), meaning that the error is proportional to h. Let y M (x m , t) ≈ y M | n m , n = 1, 2, . . . , ∞. Then, we obtain Let , and substituting (13) and (14) into (1) yields From (15), we can express the prediction value y M | n+1 m as: Since the sensor is located at x m , the approximation y M | n m can be replaced by the actual output y | n m , where the approximation y M | n+1 m is the prediction outputŷ | n+1 m . Iteratively, the recurrence formula of the P-step forward prediction output y(x, t) is given bŷ

Wavelet-Based Nonlinear Model Predictive Control
We consider (16), which is an approximation of (1). The discrete version of the objective function (2) is used to design the MPC controller. It is given by where J(∆u) is the objective function, and r | n+p m andŷ c | n+p m are, respectively, the reference trajectory and the correction prediction output defined at time t n+p and position x m . ∆u | n+i is the increment of the control signal at time n + i, u min and u max represent, respectively, the upper bound and the lower bound of the control variables, ∆u min and ∆u max represent, respectively, the upper bound and lower bound of the control increment. N 1 is the minimum cost horizon, N 2 is the prediction horizon, N c is the control horizon, λ i , i = 1, · · · , N c , are the weighting factors for the prediction error and control signal.  (20) where H is the error weight vector, which is usually chosen to be the identity vector I M×1 .
At every sample time t n , the N c -step optimal control actions ∆u * | n+1 , · · · ,∆u * | n+N c over the prediction horizon are calculated through solving the optimization problem (18). This optimization problem is, in general, much too complicated to admit an analytical solution and hence numerical solutions are to be obtained. For this, the original problem is approximated by a series of quadratic programming subproblems, and each of which is solved by using the sequential quadratic programming (SQP) algorithm [34], which is recognized as one of the most efficient algorithms for solving nonlinear optimization problems of small and moderate size. In this way, the increment of the optimal control signal is obtained. The structure of the nonlinear MPC is depicted in Figure 4.

Case Study
In this section, we consider the temperature control problem of a long duct heating system based on flow control [20], as shown in Figure 5.
This is a one-dimensional incompressible convection heat transfer process. The fluid flows from one end of the long duct to the other end. We assume that the axial direction of the long duct is the cylindrical surface of the same diameter and that the ambient temperature of the entire system, the inlet temperature and the thermal conductivity of the fluid in the long duct are constant. Furthermore, both the fluid flow and the convection heat transfer are fully developed flows, and the diffusion phenomenon is neglected. Under the above assumptions, the temperature distribution in the duct depends only on the flow velocity. This problem is divided into three categories, corresponding to the different relationships between the flow velocity and the delay time τ(t): (i). The flow velocity is high, but the delay time is short. In this case, the influence of the flow velocity on the temperature is fast; (ii). The flow velocity is slow, but the delay time is large. In this case, the influence of the flow velocity on the temperature is slow; and (iii). The flow velocity is particularly slow. In this case, the change of the temperature in response to the flow velocity is negligible. In practice, the first two treatments are similar, and the last case is rare. In actual industrial systems, the flow velocity is often affected by unknown disturbances. When the flow velocity is affected, the axial temperature of the long tube will be affected. The purpose of the controller is to quickly restore the axial temperature of the long duct back to the equilibrium state by adjusting the opening of the regulating valve. The dimensionless form of the dynamical model can be expressed as: where T(x, t) is the axis temperature of the long duct, v(t) is the fluid velocity, T t (x, t) and T x (x, t) are, respectively, the partial derivative with respect to t and x. The numerical solution is to be obtained by using FDM as described in [32]. The spatial region is evenly divided by uniform grid points {m|m = 1, 2, . . . , M}, and ∆x = 1/(M − 1). Similarly, let the time step be fixed as ∆t n , n = 1, 2, . . . , N. Then, the discretized temperature value is T | n m = T(x m , t n ) at the space location x m and the time t n . The brackets are used to distinguish between time points and indices. Let Krone number be CFL = v | n ∆t n /∆x m . Then, the discretized version of system (21) can be expressed as: The temperature value of the long duct is calculated iteratively by (22). The sensors are evenly located in a sufficiently small interval ∆x to measure the real temperature value, whereas the time interval ∆t is chosen such that Krone stability condition is satisfied, that is CFL ≤ 1. Thus, to ensure that the stability condition is satisfied, we choose ∆t ≤ ∆x/max{v | n }, where max{v | n } is the maximum flow velocity of the fluid in the long duct. The mesh is chosen to be more refine than that in [20] with the same flow velocity. Thus, the numerical solution obtained is more accurate. The length of the spatial interval ∆x is chosen to be 1/2 13 , which is small enough to minimize the inaccuracy to be an integer power of 1/2 such that it is consistent with the binary wavelet. Meanwhile, the discrete time step ∆t is chosen to be 0.0001 such that the inaccuracy is minimized and the CFL condition is satisfied. The flow velocity v(t) ∼ O(1) is selected randomly between 0.9 to 1.1. All the experiments work on IBM PC-Intel Core I5-2500K using Matlab 2016b. The actual temperature in the long duct after 10-time units is depicted in Figure 6.

The Performance of the Prediction Model
In this subsection, we will test the performance of the prediction model using the method proposed in Section 3.3. Assume that there are M = 2 j + 1 temperature sensors evenly located along the long duct. T | . The P-step prediction of the temperature value of (21) given by (17) is The choice of the parameters j, r, N of the interval quasi-Shannon interpolating wavelet will significantly affect the performance. In Table 1, we show the L ∞ -error and the mean square error (MSE) under two different sample times ∆t s = 0.1 and 0.01 for three different meshes j = 2, 3, 4 with N = 2, 3, 7. Based on the results obtained in Section 3.2, the quasi-wavelet parameters r = 2.2, 2.8, and 3.4 are chosen. Let the L ∞ error function be defined as follows: and let the MSE function be given by: From Table 1, we can see that the sample time ∆t s has a significant impact on the prediction performance. The single-step prediction performance meets the control requirements for both the cases of ∆t s = 0.1 and 0.01. On the other hand, the multi-step prediction performance fails to meet the control requirements for the case of ∆t s = 0.1. As j is increased from 2 to 3, the L ∞ error is deduced significantly. However, as j continues to increase, the decrease of the error is no longer significantly. The optimal value of the parameter r is around 2.2. The L ∞ error and MSE are the smallest for j = 2 and 3. However, the parameter r becomes insensitive to prediction performance when j = 4. The WCM, which is developed based on Haar wavelet basis in [31] and the FDM developed in [32], are selected for comparison with the method proposed in this paper. The sample time ∆t s is chosen small enough such that all errors are mainly originated from the spatial discretization. Without loss of generality, we choose ∆t s = 0.01. In Table 2, a comparison of performance between the FDM with ∆x = 1/4, 1/8, 1/16, the WCM using Haar wavelet basis with M = 4, 8, 16, and the proposed method with j = 2, 3, 4 is given. The parameters of the interval quasi-Shannon interpolating wavelet are selected as N = 2, r = 2.2 for j = 2; N = 3, r = 2.2 for j = 3; and N = 7, r = 2.8 for j = 4. The predictive temperature and the error of the model, which are obtained by the FDM with ∆x = 1/16, the WCM using Haar wavelet basis with M = 16, and the proposed method with j = 4, are depicted in Figure 7.
In Table 2, the prediction errors obtained by FDM and WCM using Haar wavelet basis decrease slowly when M increases. On the other hand, the prediction errors obtained by the proposed method decrease exponentially when j = 2 is increased to j = 3. Please note that with the increase of the integer j, the truncation error of the proposed method using the regularized Shannon's kernel decays exponentially, but the Euler forward difference method has only the first-order accuracy [35]. Therefore, we can see that the prediction performance of the model using the proposed method is better than those obtained by using the other two methods. The model based on the WCM using Haar wavelet basis with M = 4 is divergent, and the prediction performance is the worst among the three methods. From Figure 7, the error caused by using the proposed method is smaller than those induced by using the FDM and the WCM using Haar wavelet basis. The errors of the models caused by using the proposed method and the FDM are smoother than that obtained by the WCM using Haar wavelet basis. This is because Haar wavelet basis is discontinuous.

The Performance of the Proposed MPC Scheme
Based on the above prediction model, a nonlinear predictive controller is designed to deal with the temperature control problem. The controller is required to adjust the flow velocity to drive the temperature back to the desired initial steady state. Assume that the quadratic criterion to be minimized is given by (26) where J(∆v) is the objective function, and T r | n+p m andT c | n+p m are, respectively, the reference temperature value and the correction prediction temperature value defined at time t n+p and position x m . ∆v | n+i is the increment of the flow velocity at time n + i, v min and v max represent, respectively, the upper bound and the lower bound of the flow velocity, ∆v min and ∆v max represent, respectively, the upper bound and lower bound of the increment of the flow velocity. N 1 ,N 2 , N c , and λ i are defined in (18).
In this section, using the proposed method, the FDM, and the WCM using Haar wavelet basis, three sets of comparison experiments are performed. The purpose is to estimate the temperature at various spatial points along the long duct when the flow velocity v(t) is suddenly increased by 10% from 1 to 1.1 at the entrance. We suppose that there are 17 temperature sensors when the proposed method and the FDM are used, while there are 16 controlled temperature sensors when the WCM using Haar wavelet basis is used. Assume that v min = 0.9, v max = 1.1, ∆v min = −0.1 × 10 −3 , ∆v max = 0.1 × 10 −3 . The predictive horizon N 1 = 1, N 2 = 10, and the control horizon N c = 3, the weighting factors λ i = 1, i = 1, 2, 3 in the quadratic criterion expressed by (26) are chosen. The temperature response at the terminal position and the flow velocity v(t) are depicted in Figure 8, and the response curves of the temperature values at other positions x = {0.25, 0.5, 0.75, 1} are depicted in Figure 9.
From Figure 8 to Figure 9, we see that the flow velocity v(t) changes smoothly, and the temperature values at positions x = {0.25, 0.5, 0.75, 1} are driven back to the steady state position for all the three methods. The results show that all the controllers obtained can solve the temperature control problem at any measured position along the long duct. However, for the proposed method, the temperature is driven back to the steady state position at a faster rate with smaller overshoots.

Conclusions
In this paper, a first-order hyperbolic system is considered, for which an equivalent system is constructed based on the collocation technique for spatial derivatives and Euler forward difference method for the approximations of the time component. The nonlinearities are handled by the collocation method, while the accuracy of the approximate solution on the finite interval is guaranteed using the interval quasi-Shannon wavelet basis. Then, a nonlinear model predictive controller is designed through solving a nonlinear quadratic programming problem with input constraints. The FDM, the WCM using Haar wavelet basis, and the proposed method are applied to the temperature control problem of a flow-control long-duct heating system. In comparison, the multi-step prediction accuracy of the model approximation using the proposed approach is distinctly higher, and hence producing a superior control performance. From the simulation study, it is observed that the WCM has a strong potential for dealing with singular problems. For future research, it will be of interest to investigate the control of systems in the presence of local dramatic changes.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: