Multiresolution GPC-Structured Control of a Single-Loop Cold-Flow Chemical Looping Testbed

Abstract

Chemical looping is a near-zero emission process for generating power from coal. It is based on a multi-phase gas-solid flow and has extremely challenging nonlinear, multi-scale dynamics with jumps, producing large dynamic model uncertainty, which renders traditional robust control techniques, such as linear parameter varying $H_{\infty }$ design, largely inapplicable. This process complexity is addressed in the present work through the temporal and the spatiotemporal multiresolution modeling along with the corresponding model-based control laws. Namely, the nonlinear autoregressive with exogenous input model structure, nonlinear in the wavelet basis, but linear in parameters, is used to identify the dominant temporal chemical looping process dynamics. The control inputs and the wavelet model parameters are calculated by optimizing a quadratic cost function using a gradient descent method. The respective identification and tracking error convergence of the proposed self-tuning identification and control schemes, the latter using the unconstrained generalized predictive control structure, is separately ascertained through the Lyapunov stability theorem. The rate constraint on the control signal in the temporal control law is then imposed and the control topology is augmented by an additional control loop with self-tuning deadbeat controller which uses the spatiotemporal wavelet riser dynamics representation. The novelty of this work is three-fold: (1) developing the self-tuning controller design methodology that consists in embedding the real-time tunable temporal highly nonlinear, but linearly parametrizable, multiresolution system representations into the classical rate-constrained generalized predictive quadratic optimal control structure, (2) augmenting the temporal multiresolution loop by a more complex spatiotemporal multiresolution self-tuning deadbeat control loop, and (3) demonstrating the effectiveness of the proposed methodology in producing fast recursive real-time algorithms for controlling highly uncertain nonlinear multiscale processes. The latter is shown through the data from the implemented temporal and augmented spatiotemporal solutions of a difficult chemical looping cold flow tracking control problem.

1. Introduction

The current transition to clean power generation involves both the use of renewables, such as hydrokinetics [1], and cleaner coal-based techniques. The latter are projected to still supply power for the foreseeable future due to the abundance of coal in many industrialized and developing countries; however, they will be required to meet the hard caps on carbon emissions. Chemical looping (CL) is a near-zero emission coal-based technology in which multiple interacting loops of flowing, reactive, gas/solid mixtures produce energy via solid-oxygen carrier based combustion [2,3,4]. Chemical looping has remained an area of active research focused on improving its economic viability and reducing environmental footprint [5,6]. To reduce waste stream volumes and required energy, advanced optimizing control systems for the chemical looping process are required. However, the process exhibits extremely challenging nonlinear multi-scale dynamics that are hard to characterize and depend on a particular design. These features render traditional robust control techniques marginally successful in experimental trials.

The goal of the present paper is to present the development of the real-time computational control-oriented models and the corresponding model-based control design strategies found to provide the desired chemical looping tracking performance. In particular, we demonstrate a novel model-based process control methodology to control the pressure drop in the riser of a single loop chemical looping process, where the air flow rate is the controlled variable. This control approach was implemented and successfully tested on an industrial single loop cold gas/solid flow chemical looping testbed, where the previously available techniques had exhibited difficulties.

Prior to being able to control the process, it is imperative to characterize the system’s response to control inputs. Classically, this would be done by devising a physical model of the system from first principles, but this often yields limited practical utility for increasingly complex nonlinear models when viewed from the perspective of process control design. To meet this challenge, an alternative technique, identification of a model constructed on the basis of the wavelet multiresolution analysis (MRA), is used in the present work. MRA has become one of the major tools in neural networks [7,8,9,10] and nonlinear system modeling [11,12,13,14,15,16,17,18]. Wavelet-based multiresolution decomposition has been proven to constitute a universal approximator for a wide range of function spaces in terms of linear combination of scaling and wavelet functions. Wavelet approximation has no smoothness requirement on the target function, making it an appropriate candidate for identification of complex nonlinear systems with multiscale dynamics, such as those encountered in chemical looping processes. Several controller designs incorporating wavelet system representations have been proposed in the literature. Reference [18] proposed adaptive adjustment of the model resolution and the corresponding structure of the nonlinear adaptive controller. However, no optimality in controller synthesis was introduced and no testing was done on the real multiresolution system. In [19] an optimal model predictive multiresolution control law with constraints was derived. However, the controller was given as a sequence of computational steps with no clear analytical formula for controller implementation and therefore no guarantee of the acceptable real-time performance; also no adaptation was included. In [20] utilization of wavelets in generalized predictive control (GPC)l has been proposed for reduction of the computational demands on the constrained GPC, but the application was not addressing multiresolution nonlinear system modeling and was restricted to linear systems only. Thus, there has been a clearly identifiable gap in producing an optimal adaptive control law with rate constraints and guaranteed real-time performance for systems with nonlinear multiscale dynamics.

The present work fills this gap through the development of the self-tuning wavelet MRA-based topology that combines temporal and spatiotemporal loops into a single closed loop control system. The GPC structure is employed for embedding into it the identified temporal nonlinear multiscale model due to its real-timable recursive receding horizon calculation, local optimality by the virtue of being a variant of LQG [21], relatively easy incorporation of rate constraints on the control signal, tunable robustness properties (not pursued in this paper), and natural embedding of the integrator to address setpoint step changes in the chemical looping based power generation.

The paper completes the brief presentation of the results of chemical looping project at Alstom Thermal Power given in a conference publication [22] through the addition of the temporal controller derivation, and presents the previously unreleased experimental data along with the corresponding discussion, as well as the derivation and implementation of an additional spatiotemporal controller for the fast process dynamics related to the riser. The material presented is based in part on two unpublished documents, an internal technical report [23] and the PhD thesis [24] of the first author; however, the detailed controller derivation was not given in either and is presented here for the first time.

Embedding highly nonlinear, but linear in parameters, adaptable multiresolution model into GPC is a novel idea requiring nontrivial analytical effort presented in this paper. Several researchers have successfully applied GPC proposed by Clark et al. [25,26] to many control areas [27,28,29,30,31]. GPC, however, has limitations, some of which have been discussed by Grimble in [32]. Since GPC uses a linear dynamic model to make predictions of process outputs over the prediction horizon, its performance will significantly degrade when the real process has severe nonlinearities and runs in a wide range of operating conditions, as is the case for a chemical looping process. Therefore, it is imperative to incorporate a high fidelity nonlinear dynamic model into the GPC scheme. Accordingly, we embed a wavelet MRA model of the nonlinear single loop cold flow of the chemical looping process into the GPC scheme. Specifically, first, a single-input-single-output (SISO) nonlinear autoregressive exogenous model (NARX) based on wavelet MRA is identified on-line using the chemical looping process test rig. Then, a GPC-type performance index is formulated, which incorporates the MRA model, and a gradient descent (GD) algorithm is developed for tuning both the weighting parameters of the wavelet MRA model and the control sequence in the GPC scheme. Further, the Lyapunov function-based theorems are proven to separately guarantee the convergence of the wavelet MRA identified model and the stability of the proposed GPC scheme without constraint and provide a guidance on both controller and identifier performance tuning. A rate constraint is then imposed on the control signal to smooth out the CL process transients. The resulting controller is shown to yield a satisfactory closed loop performance over a broad operating range, effectively meeting the challenge of handling the chemical looping process complexity.

The resulting cold flow testbed control system was then further improved by augmenting the temporal closed loop structure described above with the additional spatiotemporal control of fast dynamics of the riser loop, which were not captured in the original low-frequency wavelet MRA system model. The response time of the 2-partial differential equations (2-PDE) riser model used for this purpose is much shorter than that of the identified NARX model, for which the sampling time is 1 s. Therefore, the control-oriented riser representation was obtained through the use of the 2-PDE riser model as follows: The model was simulated to obtain the riser impulse response, the latter was employed to approximate the faster dynamics of the system, and the result was used in a convolution to obtain a model of the transients. To simplify the calculations, the impulse response was approximated using Gaussian spatial and temporal wavelets. Simulation and experimental results verified the validity of the spatiotemporal wavelet-based control system topology augmentation.

The paper is organized as follow: Section 2.1 provides the nomenclature for the main variables and symbols used in the paper. Section 2.2 introduces the chemical looping process model. A NARX model representation and a wavelet MRA representation are given in Section 2.3. Section 2.4 provides derivation of a wavelet MRA-based GPC strategy for solving the tracking problem for a single loop cold flow system. The convergence of the prediction error of the wavelet MRA model identification algorithm and the tracking error of the proposed GPC control strategy are separately proven in Section 2.5. An input-constrained GPC scheme is presented in Section 2.6. Experimental results are discussed in Section 3. The closed loop topology augmentation with the spatiotemporal model-based control to account for the pressure drop DP47 over the riser related to the fluidizing air flows is presented in Section 4. The discussion of the results is presented in Section 5. A conclusion is provided in Section 6.

2. Materials and Methods

2.1. Nomenclature

${\varphi }_{j,k}(x)$: orthonormal basis for $V_j$; ${\psi }_{j,k}(x)$: orthonormal basis for $W_j$; $y(t)$: system output; $u(t)$: system input; $w(t)$: system noise; $e(t)$: model output estimation error; $n_y$: maximum lag in the output; $n_u$: maximum lag in the input; $\theta$: weighting parameter vector trained on-line; $g_i$: multivariable scaling or wavelet basis function of past inputs and outputs; ${\gamma }_{\theta }$: adaptation gain for the control input vector; $U$: control input vector; $\mathsf{\Delta }{u}_{\mathrm{target}}$: unconstraint control signal calculated by the predictive control law; $\mu$: design parameter; $\epsilon :$ voidage; $u_s$: solid velocity; $U_g$: superficial gas velocity; $S_{1,w}(t)$: control command calculated by wavelet adaptive GPC control; $y_r(t)$: reference signal; ${\nabla }_{\theta }{J}_1(n)$: gradient of loss function $J_1$ with respect to $\theta$; ${\nabla }_U{J}_2(n)$: gradient of the loss function $J_2$ with respect to $U$; ${\nabla }_{\theta }e(n)$: sensitivity derivative at time $n$.

2.2. Chemical Looping Process

The modeling and control methodologies proposed in this paper focus on the hybrid combustion-gasification chemical looping (CL) process initially developed by Alstom Power. Chemical looping is a two-step process which first separates oxygen (O₂) from nitrogen (N₂) in an air stream in an air reactor. The O₂ is transferred to a solid oxygen carrier. Next, the oxygen is carried by the solid oxide and is then used to gasify or combust solid fuel in a separate fuel reactor. As shown in Figure 1, a metal or calcium material (oxygen carrier) is burned in air forming a hot oxide (MeO_x or CaO_x) in the air reactor (oxidizer). The oxygen in the hot metal oxide is used to gasify coal in the fuel reactor (reducer), thereby reducing the oxide for continuous reuse in the chemical looping cycle. CL coal power technology is an entirely new, ultra clean, low cost, high efficiency coal power plant technology for the future power market. The concept promises to be the technological link from today’s steam cycle power plants to tomorrow’s clean coal power plants, capable of high efficiency and CO₂ capture.

The CL process with its multi-phase flows and complicated chemical reactions is characterized by process nonlinearities and time delays due to mass transport and chemical reactions. The specific operational characteristics are new and are still being studied. Hence, there is a need for further investigation and the potential for advanced control solutions. In this paper, we have focused on developing a control-oriented model for a single loop cold gas/solid flow test rig which omits all chemical reactions and interactions with other loops.

The block diagram of a single loop cold flow CL process is shown in Figure 2. It consists of a lower level pipeline, a riser pipeline, an upper level horizontal pipeline, a cyclone, a dip leg, seal pot control valves (SPCV), and a solid return leg. The lower level pipeline accepts air flow and solids returned from both seal pot control valves and solids which are added manually. In the riser the air-solid mixture (two-phase) flows upwards, turns into the horizontal pipeline, and then enters the cyclone. The cyclone separates the solid particles from the air. The separated solids then drop into the dip leg and enter the SPCV. The SPCV splits the solids between the return leg in its own loop and the return leg in another loop. The SPCV also maintains a pressure control boundary.

In our test rig, the manipulated variables (MV) include $S1,S2$—two fluidizing air flow rates into the SPCV, which change pressures in the SPCV and the flow conditions upstream and downstream of the SPCV. The controlled variable (CV) of interest is $\mathrm{DP}47$, which stands for the pressure drop measured across the riser—a substantive indicator of solid/gas flow transport stability along the whole loop. The performance of the test rig implementing the controller to track a reference command was evaluated both under step-changes and cycling operation.

A wavelet MRA modelling and its embedding into a GPC-based predictive controller are described in the next two subsections.

2.3. Wavelet MRA Model Structure

Wavelet multiresolution analysis [14] is a function approximation tool representing function details at different scales of resolution in both the time and the frequency domains in terms of shifted and dilated scaling and wavelet functions. In general, MRA consists of a sequence of successive approximation closed subspaces $V_j\in L_2(R)$, $j\in Z$ satisfying:

$$\cdots V_{-1}\subset V_0\subset V_1\cdots ,$$

(1)

with the following properties:

$${\cup }_{j\in Z}{V}_j \mathrm{is} \mathrm{dense} \mathrm{in} L_2(R); {\cap }_{j\in Z}{V}_j=\left\{0\right\},$$

(2)

$$f(x)\in V_j\iff f(2x)\in V_{j+1}$$

(3)

$$f(x)\in V_j\iff f(x-2^{-j}k)\in V_j, k\in Z,$$

(4)

$$V_j=\mathrm{span}\{{\varphi }_{j,k},k\in Z\},$$

(5)

where $Z$ is the set of all integers, ${\varphi }_{j,k}(x)=2^{j/2}\varphi (2^{j}x-k)$ is an orthonormal basis for $V_j$ and $L_2(R)$ is the space of square integrable functions of scalar real variable.

If we define $W_j$ to be the orthogonal complement of $V_j$ in $V_{j+1}$, then:

$$V_{j+1}=V_j\oplus W_j, V_j\perp W_j,$$

(6)

$$W_j=\mathrm{span}\{{\psi }_{j,k},k\in Z\},$$

(7)

where ${\psi }_{j,k}(x)=2^{j/2}\psi (2^{j}x-k)$ is an orthonormal basis for $W_j$. It follows from Equations (1) and (6) that, any $V_j$ can be written for any $l<j$ as:

$$V_{j+1}=V_l\oplus W_l\oplus W_{l+1}\oplus W_{l+2}\oplus \cdots \oplus W_j,$$

(8)

where all the subspaces are orthogonal. Then this implies that:

$$L_2(R)={\oplus }_{j\in Z}{W}_j$$

(9)

The functions ${\varphi }_{j,k}$ and ${\psi }_{j,k}$ will be referred to as scaling and wavelet functions respectively. According to Equations (8) and (9), any $f(x)\in L_2(R)$ can be represented as:

$$f(x)={\displaystyle \sum {}_n(f,{\varphi }_{J,n})}{\varphi }_{J,n}+{\displaystyle \sum {}_{j\ge J,n}(f,{\psi }_{j,n}){\psi }_{j,n}}$$

(10)

The approximation starts from some lower resolution level $J$ and can be truncated at certain higher resolution level $N$ when:

$$\parallel f(x)-[{\displaystyle \sum _n(f,{\varphi }_{J,n})}{\varphi }_{J,n}(x)+{\displaystyle \sum _{j=J}^N{\displaystyle \sum _n(f,{\psi }_{j,n})}}{\psi }_{j,n}(x)]\parallel <\epsilon ,$$

(11)

for any predefined small error $\epsilon >0$.

Multivariable wavelet bases can be constructed from the tensor product of a radial basis function of a one-dimensional wavelet as described for images in [33]. Because wavelet MRA can approximate any finite energy nonlinear function to any desired accuracy level, in this paper, the wavelet MRA will be used to build the nonlinear empirical model for a single loop cold flow CL process, as shown in the next subsection.

The NARX Model Structure

Many systems in a variety of applications contain nonlinearities which render linear model incapable of capturing the complex dynamic system behavior. Therefore, it is of interest to develop for these applications sufficiently accurate nonlinear dynamical models. An NARX model [34] is a well-established input/output representation for nonlinear system identification. Under some mild assumptions, a discrete-time stochastic nonlinear SISO system can be expressed as:

$$y(t)=f(y(t-1),\cdots ,y(t-n_y),u(t-1),\cdots ,u(t-n_u))+w(t),$$

(12)

where $y(t)$, $u(t)$, $w(t)$ are the system output, input, and noise, and $t$ is discrete time, respectively, $n_y$ and $n_u$ are the maximum lags in the output and input, $w(t)$ is assumed to be a zero mean, independent, and bounded noise variable, and $f(\cdot )$ is some nonlinear function. Unless some prior knowledge of the system dynamics is available, most methods use nonparametric regression to estimate the nonlinear function $f$ from the data. In our case, $f$ is implemented as a linear expansion in terms of the scaling and wavelet functions of regressors $g_i$ such that

$$f={{\displaystyle \sum }}_{i=1}^m{\theta }_i{g}_i$$

(13)

minimizes a pre-specified approximation adequacy criterion, where $\theta =\{{\theta }_i\}$ is a parameter vector trained on-line, $g_i\in \{{\varphi }_{j,k},{\psi }_{j,k}\}$ is a multivariable scaling or wavelet basis function of past inputs and outputs, and $m$ is the number of basis functions to meet some given modeling accuracy requirement.

In the next subsection, the NARX model structure introduced above is embedded into the parameter adaptation law and the GPC performance criterion, and the self-tuning MRA-based control law is derived.

2.4. Wavelet MRA-Based GPC Scheme

The basic methodology of GPC is to calculate the current control actions on-line at each sampling instant in order to solve a finite horizon, open-loop, optimal control problem where the first control in the optimal control sequence is applied to the plant. In this section, we present both the online wavelet MRA system identification algorithm and the GPC based predictive control strategy for the stable tracking problem of a single loop CL system. To clearly illustrate the idea of the proposed control scheme, we derive the algorithm for a SISO nonlinear dynamic system. The extension to a multi-input- multi-output setting is straightforward.

Referring to Figure 2 and its description, let $\mathrm{DP}47$ be the actual system output $y$ and $S1$ be the control input $u$, while $S2$ is set to a constant value. Let $\widehat{y}$ denote the approximated system output. Then, the identified wavelet MRA based model is defined as follows:

$$\widehat{y}(t)=f(y(t-1),\cdots ,y(t-n_y),u(t-1),\cdots ,u(t-n_u)),$$

(14)

where f is defined in Equation (13). Then, the error between the real plant output $y$ and the estimated output $\widehat{y}$ is defined as:

$$e(n)=y(n)-\widehat{y}(n).$$

(15)

The weighting parameters $\theta$ in Equation (13) are trained online to minimize the loss function defined as:

$$J_1(n)=\frac{1}{2}{e}^2(n),$$

(16)

where $n$ indicates discrete time. To make $J_1$ small, we employ a parameter adaptation law in the form of a gradient descent (GD) algorithm, which adjusts the weighting gains $\theta$ to keep the gradient of $J_1$ negative, that is:

$$\theta (n+1)=\theta (n)-{\gamma }_{\theta }{\nabla }_{\theta }{J}_1(n)=\theta (n)-{\gamma }_{\theta }e(n){\nabla }_{\theta }e(n),$$

(17)

where ${\gamma }_{\theta }$ is the adaptation gain, ${\nabla }_{\theta }{J}_1(n)$ is the gradient of $J_1$ with respect to $\theta$ at discrete time $n$, and ${\nabla }_{\theta }e(n)$ is the so-called sensitivity derivative at time $n$ indicating how the error is influenced by the weighting parameters $\theta$. From Equations (13)–(15), the sensitivity derivative ${\nabla }_{\theta }e$ can be derived as follows:

$${\nabla }_{\theta }e=-{\nabla }_{\theta }\widehat{y}=-{\nabla }_{\theta }f=-g\Rightarrow \theta (n+1)=\theta (n)+{\gamma }_{\theta }e(n)g(n).$$

(18)

Suppose the future set-point signals $y_m(n+k),k=1,2,\cdots$ are available. In the context of GPC, define another loss function as follows:

$$J_2=\frac{1}{2}\{{{\displaystyle \sum }}_{k=N_1}^{N_2}{(y_m(n+k)-\widehat{y}(n+k))}^2+{{\displaystyle \sum }}_{k=1}^{N_u}{\rho }_k\mathsf{\Delta }u{(n+k-1)}^2\},$$

(19)

where $N_1$ and $N_2$ are the minimum and the maximum output prediction horizons, respectively, $N_u$ is the control horizon, $\mathsf{\Delta }$ is the difference operator, $\mathsf{\Delta }u(n)=u(n)-u(n-1)$, and ${\rho }_k$ is the $k$-th control weighting factor. Assuming $N_1=1$, $N_2=L=N_p$, and identical control weighing factor ${\rho }_k=\rho ,$ Equation (19) can be rewritten in the vector form as:

$$J_2=\frac{1}{2}\left\{{\parallel Y_m(n+1)-\widehat{Y}(n+1)\parallel }^2+\rho {\parallel \mathsf{\Delta }U(n)\parallel }^2\right\},$$

(20)

where:

$$Y_m(n+1)={[y_m(n+1),y_m(n+2),\cdots ,y_m(n+L)]}^T,$$

$$\widehat{Y}(n+1)={[\widehat{y}(n+1),\widehat{y}(n+2),\cdots ,\widehat{y}(n+L)]}^T,$$

$$U(n)={[u(n),u(n+1),\cdots ,u(n+N_u-1)]}^T,$$

$$\mathsf{\Delta }U(n)={[\mathsf{\Delta }u(n),\mathsf{\Delta }u(n+1),\cdots ,\mathsf{\Delta }u(n+N_u-1)]}^T,$$

and $\parallel \cdot \parallel$ is the $L_2$ norm of the $n$-dimensional real vectors.

The loss function $J_2$ is now minimized to drive the system output $y$ to the reference signal $y_m$ given that the wavelet MRA identifier accurately approximates the real process dynamics on-line. At each sampling instant, an optimal control sequence is calculated using future predicted output values of the identified model, but only the first one is applied to the system. To minimize $J_2$, the GD method is implemented again to recursively calculate the $N_u$-dimensional control increment sequence $\mathsf{\Delta }U$ as follows:

$$\mathsf{\Delta }U(n)=-{\gamma }_u{\nabla }_U{J}_2(n),$$

(21)

where ${\gamma }_u$ is the adaptation gain for the control input vector $U$. Noting that for any vector y(x), ${\nabla }_x{\parallel y\parallel }^2=2({\nabla }_{x}y)y$, from Equations (20) and (21), the gradient of the loss function $J_2$ with respect to $U$ can be obtained as:

$$\begin{array}{cc}\hfill {\nabla }_U{J}_2(n)& =\frac{1}{2}{\nabla }_U\left\{{\parallel Y_m(n+1)-\widehat{Y}(n+1)\parallel }^2+\rho {\parallel \mathsf{\Delta }U(n)\parallel }^2\right\}\hfill \\ & ={\nabla }_U\left\{Y_m(n+1)-\widehat{Y}(n+1)\right\}\left\{Y_m(n+1)-\widehat{Y}(n+1)\right\}+\rho {\nabla }_U\left\{\mathsf{\Delta }U(n)\right\}\left\{\mathsf{\Delta }U(n)\right\}.\hfill \end{array}$$

The first part of the expression above is evaluated as follows. First, we note that since $y_m(n)$ are the reference signals, which are preset constants, we have $\partial y_m(n+1)/\partial u(n)=0$. Then, since $\widehat{y}(n+k)={\theta }^{T}u$, $\widehat{y}$ depends only on the past u, we have:

$$\frac{\partial y(n+k)}{\partial u(n+l)}=\{\begin{array}{c}\frac{\partial \widehat{y}(n+k)}{\partial u(n+l)}, \mathrm{when} k>l,\\ 0, \mathrm{when} k\le l.\end{array}$$

This yields:

$${\nabla }_U\left\{Y_m(n+1)-\widehat{Y}(n+1)\right\}=$$

$${\left[\begin{array}{cccc}\frac{\partial }{\partial u(n)}\left\{y_m(n+1)-\widehat{y}(n+1)\right\}& \frac{\partial }{\partial u(n)}\left\{y_m(n+2)-\widehat{y}(n+2)\right\}& \cdots & \frac{\partial }{\partial u(n)}\left\{y_m(n+L)-\widehat{y}(n+L)\right\}\\ \frac{\partial }{\partial u(n+1)}\left\{y_m(n+1)-\widehat{y}(n+1)\right\}& \frac{\partial }{\partial u(n+1)}\left\{y_m(n+2)-\widehat{y}(n+2)\right\}& \cdots & \frac{\partial }{\partial u(n+1)}\left\{y_m(n+L)-\widehat{y}(n+L)\right\}\\ ⋮& ⋮& \ddots & ⋮\\ \frac{\partial }{\partial u(n+N_u-1)}\left\{y_m(n+1)-\widehat{y}(n+1)\right\}& \frac{\partial }{\partial u(n+N_u-1)}\left\{y_m(n+2)-\widehat{y}(n+2)\right\}& \cdots & \frac{\partial }{\partial u(n+N_u-1)}\left\{y_m(n+L)-\widehat{y}(n+L)\right\}\end{array}\right]}_{N_u\times N_p}$$

$$={\left[\begin{array}{cccccc}-\frac{\partial \widehat{y}(n+1)}{\partial u(n)}& -\frac{\partial \widehat{y}(n+2)}{\partial u(n)}& \begin{array}{ccc}& & \cdots \end{array}& -\frac{\partial \widehat{y}(n+N_u)}{\partial u(n)}& \begin{array}{cc}\cdots & \end{array}& -\frac{\partial \widehat{y}(n+L)}{\partial u(n)}\\ 0& -\frac{\partial \widehat{y}(n+2)}{\partial u(n+1)}& \begin{array}{ccc}& & \cdots \end{array}& -\frac{\partial \widehat{y}(n+N_u)}{\partial u(n+1)}& \begin{array}{cc}\cdots & \end{array}& -\frac{\partial \widehat{y}(n+L)}{\partial u(n+1)}\\ & & & \begin{array}{ccc}⋮& & \end{array}& & \\ \begin{array}{ccccc}& & 0& & \cdots \end{array}& 0& & \begin{array}{cccc}-\frac{\partial \widehat{y}(n+N_u)}{\partial u(n+N_u-1)}& & & \cdots \end{array}& & -\frac{\partial \widehat{y}(n+L)}{\partial u(n+N_u-1)}\end{array}\right]}_{N_u\times N_p}$$

The second part of ${\nabla }_U{J}_2(n)$ is evaluated by taking into account the relation $\mathsf{\Delta }u(n)=u(n)-u(n-1),$ so that $\partial \mathsf{\Delta }u(n)/\partial u(n)=1$ and $\partial \mathsf{\Delta }u(n)/\partial u(n-1)=-1$. The latter yields:

$${\nabla }_U\left\{\mathsf{\Delta }U(n)\right\}={\left[\begin{array}{ccccc}\frac{\partial \mathsf{\Delta }u(n)}{\partial u(n)}& \frac{\partial \mathsf{\Delta }u(n+1)}{\partial u(n)}& \frac{\partial \mathsf{\Delta }u(n+2)}{\partial u(n)}& \cdots & \frac{\partial \mathsf{\Delta }u(n+N_u-1)}{\partial u(n)}\\ \frac{\partial \mathsf{\Delta }u(n)}{\partial u(n+1)}& \frac{\partial \mathsf{\Delta }u(n+1)}{\partial u(n+1)}& \frac{\partial \mathsf{\Delta }u(n+2)}{\partial u(n+1)}& \cdots & \frac{\partial \mathsf{\Delta }u(n+N_u-1)}{\partial u(n+1)}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ \frac{\partial \mathsf{\Delta }u(n)}{\partial u(n+N_u-1)}& \cdots & \frac{\partial \mathsf{\Delta }u(n+N_u-3)}{\partial u(n+N_u-1)}& \frac{\partial \mathsf{\Delta }u(n+N_u-2)}{\partial u(n+N_u-1)}& \frac{\partial \mathsf{\Delta }u(n+N_u-1)}{\partial u(n+N_u-1)}\end{array}\right]}_{N_u\times N_u}$$

(22)

$$={\left[\begin{array}{ccccc}1& -1& 0& \cdots & 0\\ 0& 1& -1& \cdots & 0\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ 0& \cdots & 0& 0& 1\end{array}\right]}_{N_u\times N_u}.$$

Combining the above expressions into:

$${\nabla }_U{J}_2(n)=-G\left(Y_m(n+1)-\widehat{Y}(n+1)\right)+\rho H\mathsf{\Delta }U(n),$$

and substituting into Equation (21) as:

$$\mathsf{\Delta }U(n)=-{\gamma }_u{\nabla }_U{J}_2(n)={\gamma }_{u}G\left(Y_m(n+1)-\widehat{Y}(n+1)\right)-{\gamma }_u\rho H\mathsf{\Delta }U(n),$$

(23)

where:

$$G={\left[\begin{array}{cccccc}\frac{\partial \widehat{y}(n+1)}{\partial u(n)}& \frac{\partial \widehat{y}(n+2)}{\partial u(n)}& \begin{array}{ccc}& & \cdots \end{array}& \frac{\partial \widehat{y}(n+N_u)}{\partial u(n)}& \begin{array}{cc}\cdots & \end{array}& \frac{\partial \widehat{y}(n+L)}{\partial u(n)}\\ 0& \frac{\partial \widehat{y}(n+2)}{\partial u(n+1)}& \begin{array}{ccc}& & \cdots \end{array}& \frac{\partial \widehat{y}(n+N_u)}{\partial u(n+1)}& \begin{array}{cc}\cdots & \end{array}& \frac{\partial \widehat{y}(n+L)}{\partial u(n+1)}\\ & & & \begin{array}{ccc}⋮& & \end{array}& & \\ \begin{array}{ccccc}& & 0& & \cdots \end{array}& 0& & \begin{array}{cccc}\frac{\partial \widehat{y}(n+N_u)}{\partial u(n+N_u-1)}& & & \cdots \end{array}& & \frac{\partial \widehat{y}(n+L)}{\partial u(n+N_u-1)}\end{array}\right]}_{N_u\times N_p}$$

(24)

and:

$$H={\left[\begin{array}{ccccc}1& -1& 0& \cdots & 0\\ 0& 1& -1& \cdots & 0\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ 0& \cdots & 0& 0& 1\end{array}\right]}_{N_u\times N_u}$$

(25)

yields the control law of the form:

$$\mathsf{\Delta }U(n)={(I+{\gamma }_u\rho H)}^{-1}{\gamma }_{u}G\left(Y_m(n+1)-\widehat{Y}(n+1)\right),$$

(26)

where $I$ is the $N_u\times N_u$ identity matrix. $G$ can be computed from the chosen wavelet MRA model structure. The proposed wavelet MRA model-based GPC control schematic is shown in Figure 3. As a result, the tracking problem for a single loop cold flow system can be solved by the wavelet MRA-based GPC control strategy using the convergence tuning guidelines developed in the next section.

2.5. Convergence and Stability

In this section, we show the output error convergence of the wavelet MRA model identification algorithm and the tracking error convergence of the proposed GPC-based control strategy. These proofs serve to show that the system identification scheme will converge to the true system model of the preselected resolution, while the predictive control scheme will provide tracking of the desired output by the system output. The adaptive identification and control laws have one parameter each in the form of the adaptation gains chosen by the user. It has been shown [35] that adaptation gains are crucial to the stability and performance of an adaptive control system. Therefore, we have provided analytic guidelines for selecting these gains. The validity of such separate convergence analysis is certainly limited under significant coupling between identification and control (e.g., for aggressively chosen gains), and the coupled analysis will be reported elsewhere. However, the results are well supported by the actual implementation and testing.

2.5.1. Convergence of Wavelet MRA Identifier

Define a discrete-type Lyapunov function as:

$$V_1(n)=\frac{1}{2}{e}^2(n),$$

(27)

where $e(n)$ defined in Equation (15) represents the output modeling error. Then, the increment of the Lyapunov function is given by:

$$\mathsf{\Delta }{V}_1(n)=V_1(n+1)-V_1(n)=\frac{1}{2}\left(e^2(n+1)-e^2(n)\right).$$

(28)

The error difference can be represented using the Jacobian matrix by:

$$\mathsf{\Delta }e(n)=e(n+1)-e(n)=\left[\frac{\partial e(n)}{\partial \theta (n)}\right]\mathsf{\Delta }\theta (n)$$

(29)

where $\mathsf{\Delta }\theta (n)={\{\mathsf{\Delta }{\theta }_i(n)\}}_{i=1}^m$ represents a change in the arbitrary component of the weighting gain vector $\theta$. From Equation (18), $\mathsf{\Delta }{\theta }_i(n)$ can be obtained by:

$$\mathsf{\Delta }\theta (n)={\gamma }_{\theta }e(n)g(n),$$

(30)

$$\left[\frac{\partial e(n)}{\partial \theta (n)}\right]=-\left[\frac{\partial \widehat{y}(n)}{\partial \theta (n)}\right]=-g^T(n),$$

(31)

where $g(n)=\hspace{0.17em}{\{g_i(n)\}}_{i=1}^m$.

Theorem 1.

Let${\gamma }_{\theta }$be the adaptation gain for the weights of the wavelet MRA identified model and$g_{\max }$be defined as$g_{\max }:={\max }_n\parallel g(n)\parallel$, where$g$is the wavelet MRA basis function,$n$is the discrete time index, and$\parallel \cdot \parallel$is the$L_2$norm of a real vector. Then convergence is guaranteed if${\gamma }_{\theta }$is chosen as:

$$0<{\gamma }_{\theta }<\frac{2}{g_{\max }^2}$$

(32)

Proof. 

From Equations (28)–(31), $\mathsf{\Delta }{V}_1(n)$ can be represented as:

$$\mathsf{\Delta }{V}_1(n)=\mathsf{\Delta }e(n)\left[e(n)+\frac{1}{2}\mathsf{\Delta }e(n)\right]=\left[\frac{\partial e(n)}{\partial \theta }\right]{\gamma }_{\theta }e(n)g(n)\times \left\{e(n)+\frac{1}{2}\left[\frac{\partial e(n)}{\partial \theta }\right]{\gamma }_{\theta }e(n)g(n)\right\}$$

$$=-{\gamma }_{\theta }{e}^2(n){\parallel g(n)\parallel }^2+\frac{1}{2}{\gamma }_{\theta }^2{e}^2(n){\parallel g(n)\parallel }^4=-\lambda e^2(n),$$

(33)

where:

$$\lambda =\frac{1}{2}{\gamma }_{\theta }{\parallel g(n)\parallel }^2(2-{\gamma }_{\theta }{\parallel g(n)\parallel }^2)\ge \frac{1}{2}{\gamma }_{\theta }{\parallel g(n)\parallel }^2(2-{\gamma }_{\theta }{g}_{\max }^2)>0.$$

(34)

From Equation (32) we obtain $V_1(n)\ge 0$ and $\mathsf{\Delta }{V}_1(n)<0$, then the convergence of the weighting parameters of the identified wavelet MRA model is guaranteed. ☐

2.5.2. Stability Analysis of Wavelet MRA-Based GPC

Define a second discrete Lyapunov function as:

$$V_2(n)=\frac{1}{2}{\parallel E(n+1)\parallel }^2,$$

(35)

where $E(n+1)=Y_m(n+1)-\widehat{Y}(n+1)$. Then the change of the Lyapunov function is obtained as:

$$\mathsf{\Delta }{V}_2(n)=V_2(n+1)-V_2(n)=\frac{1}{2}({\parallel E(n+2)\parallel }^2-{\parallel E(n+1)\parallel }^2).$$

(36)

Similarly to Equation (29), the error difference can be represented using the Jacobian matrix by:

$$\mathsf{\Delta }E(n+1)=E(n+2)-E(n+1)=\left[\frac{\partial E(n+1)}{\partial U(n)}\right]\mathsf{\Delta }U(n),$$

(37)

where $\mathsf{\Delta }U(n)$ is defined in Equation (26) and $\frac{\partial E(n+1)}{\partial U(n)}=-G^T$. Then Equation (37) can be expressed as:

$$\mathsf{\Delta }E(n+1)=-G^T{(I+{\gamma }_u\rho H)}^{-1}{\gamma }_{u}GE(n+1).$$

(38)

Theorem 2.

Let${\gamma }_u$be the adaptation gain for the GPC control input sequence. Assume a control weighting factor$\rho >0$. Then the stable tracking convergence of the wavelet MRA based GPC control system is guaranteed if:

$$0<{\gamma }_u<\frac{2}{{\lambda }_{\max }(G{G}^T)},$$

(39)

where${\lambda }_{\max }(\cdot )$is the maximum eigenvalue of the matrix.

Proof. 

From Equations (36)–(38), $\mathsf{\Delta }{V}_2(n)$ can be represented as:

$$\begin{array}{c}\mathsf{\Delta }{V}_2(n)=\frac{1}{2}\left[{(E(n+1)+\mathsf{\Delta }E(n+1))}^T(E(n+1)+\mathsf{\Delta }E(n+1))-E{(n+1)}^{T}E(n+1)\right]\hfill \\ =\mathsf{\Delta }{E}^T(n+1)\left[E(n+1)+\frac{1}{2}\mathsf{\Delta }E(n+1)\right]\hfill \\ =-{(GE)}^T{\gamma }_u{({(I+{\gamma }_u\rho H)}^{-1})}^T\left[I-\frac{1}{2}G{G}^T{(I+{\gamma }_u\rho H)}^{-1}{\gamma }_u\right]GE\hfill \\ =-{(GE)}^T{R}_1{R}_{2}GE,\hfill \end{array}$$

(40)

where:

$$R_1={\gamma }_u{({(I+{\gamma }_u\rho H)}^{-1})}^T,$$

(41)

$$R_2=I-\frac{1}{2}{\gamma }_{u}G{G}^T{(I+{\gamma }_u\rho H)}^{-1}.$$

(42)

 ☐

If $R_1$ and $R_2$ are both positive definite matrices, it follows that $\mathsf{\Delta }{V}_2(n)<0$. Together with $V_2(n)>0$, the stable tracking of the reference signals is guaranteed.

From Equation (25) it can be shown that the eigenvalues of $H$ are ${\lambda }_H={\{1,\cdots ,1\}}_{N_u\times 1}$. Then the eigenvalues of $R_1$ can be derived as:

$${\lambda }_{R_1}={\{{\gamma }_u{(1+{\gamma }_u\rho )}^{-1},\cdots ,{\gamma }_u{(1+{\gamma }_u\rho )}^{-1}\}}_{N_u\times 1}$$

(43)

Hence, all eigenvalues of $R_1$ are positive if ${\gamma }_u>0$. It follows that $R_1>0$.

If $0<{\gamma }_u<\frac{2}{{\lambda }_{\max }(G{G}^T)}$, then:

$$I-\frac{1}{2}{\gamma }_{u}G{G}^T>0.$$

(44)

From Equation (25) we have:

$${\gamma }_u\rho H>0.$$

(45)

Then from Equations (44) and (45)

$$I-\frac{1}{2}{\gamma }_{u}G{G}^T+{\gamma }_u\rho H>0.$$

(46)

Similarly to the way it was done for Equation (43), we can prove that $I+{\gamma }_u\rho H>0$. Then Equation (46) can be rewritten as:

$$(I+{\gamma }_u\rho H)\left(I-\frac{1}{2}{\gamma }_{u}G{G}^T{(I+{\gamma }_u\rho H)}^{-1}\right)>0.$$

(47)

Since $I+{\gamma }_u\rho H>0$, $\left(I-\frac{1}{2}{\gamma }_{u}G{G}^T{(I+{\gamma }_u\rho H)}^{-1}\right)>0$ follows. Now we have $V_2(n)\ge 0$ and $\mathsf{\Delta }{V}_2(n)<0$. With this, the convergence of the prediction error of the wavelet MRA model identification algorithm and the tracking error of the proposed GPC control strategy have been separately proven.

2.6. Wavelet MRA GPC with Input Constraints

The stability analysis in Section 2.5 does not account for constraints. In practice, all process inputs are subject to certain constraints due to the actuation limits. In [36], two types of constraints are considered in the GPC design procedure, namely the rate and the magnitude limits on the input control signal, given, respectively, by:

$$\mathsf{\Delta }{u}_{\min }\le u(n+k)-u(n+k-1)\le \mathsf{\Delta }{u}_{\max },$$

(48)

$$u_{\min }\le u(n+k)\le u_{\max },$$

(49)

where $0\le k\le N_u-1$. When constraints are included, the stability properties obtained above must be reanalyzed. The stability analysis for constrained wavelet MRA–GPC architecture will be addressed elsewhere. Taking into account the CL process actuator constraints, the control input $u$ is subject to an input magnitude constraint saturation:

$$u^{*}(n)=\mathrm{sat}[u(n)]=\{\begin{array}{c}{u}_{\min } \mathrm{if} u(n)<u_{\min }\\ u_{\max } \mathrm{if} u(n)>u_{\max }\\ u(n) \mathrm{otherwise}\end{array}.$$

(50)

In the experiments, the latter constraints were seldomly attained, whereas control rate constraints of the form of Equation (48) had to be introduced to achieve good experimental results, as presented in the next section.

3. Results

3.1. MRA Temporal Modeling of the Chemical Looping Process Testbed

This section describes implementation of the proposed wavelet MRA model-based GPC scheme on the single loop gas/solid cold flow CL process testbed developed at Alstom Power Inc. to carry out experiments without consideration of the oxidation reaction. The experimental facility is shown in Figure 4.

The controllers were developed in MATLAB/SIMULINK, compiled in C and run on the proprietary ASTOM processing platform. The software used for wavelet identification was MATLAB Wavenet. The system output $y$ was selected to be the riser pressure drop DP47 (inch H₂O). Fluidizing air flow $S1$ (standard cubic feet per hour, scfh), was used as the single control input u, while the other air flow $S2$ (scfh) was set to a constant value of about 20 scfh. The characterization of the complex dynamic behavior, to be obtained through the identification procedure, was chosen as a SISO NARX wavelet multiresolution network model of the form:

$$\widehat{y}(t)=f(y(t-1),\cdots ,y(t-n_y),u(t-1),\cdots ,u(t-n_u))={\displaystyle \sum _{i=1}^m{\theta }_i{g}_i},$$

(51)

where $f$ is the unknown nonlinear mapping to be identified, $u(t)$ and $y(t)$ are the sampled input and output sequences, $n_y$ and $n_u$ are the maximum lags in the output and the input to be determined, respectively; $\theta =\{{\theta }_i\}$ is the parameter vector trained on-line, $g_i\in \{{\varphi }_{j,k},{\psi }_{j,k}\}$ is a multivariable scaling or wavelet basis function of past inputs and outputs, and $m$ is the number of required basis functions to meet satisfactory modeling accuracy requirements.

First, several offline experimental tests were performed to understand the process better and to leverage the test results in tuning the identification structure and the model parameters. The input signal $S_1$ was generated in the form of a pseudo random binary signal (PRBS) changed about a nominal value, and the pressure drop across the riser $\mathrm{DP}47$ was measured as the output. All the sequences used in the experiment were generated by MATLAB commands. The experimental results of the PRBS test are shown in Figure 5 and Figure 6 below, where the sampling period is 1 s.

The data set consists of 3961 input and output samples. The NARX multiresolution wavelet MRA model was used to approximate the nonlinear relationship between $S1$ and $\mathrm{DP}47$ based on the experimental data. The regressor set was specified as:

$$y(t-1),y(t-2),u(t-1),u(t-2),\cdots ,u(t-8).$$

(52)

Hence $n_y=2, n_u=8$. For MRA model we chose radial Marr scaling and wavelet functions [12]:

$$\varphi (x)=\exp (-0.5{\parallel x\parallel }^2), \psi (x)=(\dim (x)-{\parallel x\parallel }^2)\exp (-0.5{\parallel x\parallel }^2).$$

(53)

The initial coarse layer index $J$ was chosen to be 3, with the number of basis functions doubling when resolution was incremented by 1 starting with 10. The final resolution adopted was $K=6$. Figure 7 below shows how the model predicted output compares with the experimental results.

The one-step-ahead predicted output and the test data set are shown in Figure 8. From these figures it is seen that the NARX wavelet MRA model obtained predicted the system outputs well. The model was found to be sufficiently accurate and no finer resolution levels were needed to be added to the model structure.

3.2. MRA–GPC Scheme Implementation on a Chemical Looping Process Testbed

In order to make the wavelet MRA GPC applicable to the CL process, the control input $u$ was subjected to rate constraints of the form:

$$\mathsf{\Delta }u=\mathsf{\Delta }{u}_{\mathrm{target}}\times \exp (1-\mu \parallel \mathsf{\Delta }{u}_{\mathrm{target}}\parallel ),$$

(54)

where $\mathsf{\Delta }{u}_{\mathrm{target}}$ is the unconstraint control signal calculated by the predictive control law and $\mu >0$ is a design parameter to adjust the rate of the control signal. The effectiveness of such input constrained wavelet predictive controller on the CL process is demonstrated next through experimental results.

The control objective of the GPC design is to ensure that the output of the system $y$ asymptotically tracks the reference signal $y_m$. The cost function to be minimized is defined in (19). The design parameters for GPC configuration were chosen as $N_1=1, N_2=10, N_u=8, \rho =1$. The adaptation gains derived from Theorems 1 and 2, were chosen as ${\gamma }_{\theta }=0.01, {\gamma }_u=0.1$. The system was initially stable around level of $y_0=13$ inch H₂O. Two setpoint step change experiments were performed consecutively. After 5 min, the setpoint was first increased to $16$ inch H₂O and stayed at the latter value for about 7 min. Then it went back to the original level of $13$ inch H₂O. The air flow $S2$ was set to a constant value of about 20 scfh. The tracking response of the system output and the corresponding control efforts are shown in Figure 9 and Figure 10, respectively. It can be seen from these figures that the proposed wavelet MRA-based GPC method effectively tracks the setpoint changes for a single loop CL process.

As can be seen in Figure 9, starting at a setpoint of 13 inH₂O for the duration of five minutes, responds to a setpoint change at 14:20 within roughly two minutes. While there is an initial overshoot, this overshoot has a magnitude of one inH₂O, and is quickly subdued; a similar phenomenon is seen at 14:27, when the setpoint is restored to 13 inH₂O. This points to an adequate response to step-changes in $\mathrm{DP}47$ setpoint, which the process control methodology is capable of handling due to a sufficiently accurate wavelet MRA model, and effective GPC tuning; the latter can be seen in Figure 10, where controller rates are limited to within feasible bounds, and control efforts are limited, making for a subdued control input history.

In the second test, reference signal was set to be a sinusoidal of the form $y_m(t)=13+2\sin (2\pi \times 0.01\times t)$, while $S2$ (scfh) was still held at a constant value of about 20. The tracking response of system output and the corresponding control effort, shown in Figure 11 and Figure 12, respectively, demonstrate that the controller satisfies the tracking performance requirement, with a time delay between the control signal and system output potentially addressed through prediction adjustment.

In particular, a relatively timely response to reference signal changes is clearly seen in Figure 11; the controller effectuates the reference signal in about 50 s, with little overshoot, as was seen in the previous setpoint step-change experiment. However, besides the phase difference between the reference and true output, large values of undershoot are seen. This can be attributed to the presence of overshoot; as GPC reacts to samples 10 s ahead of time, any overshoot is met with an overaggressive response to lower it (see Figure 12), resulting in excessive undershoot. This issue may be addressed by increasing the control weighing factor to penalize excessive actuation.

The next section presents the derivation and implementation of the spatiotemporal control law for the fast riser dynamics to augment the temporal controller described above and tighten the closed loop tracking performance.

4. Spatiotemporal Wavelet Decomposition

Since the empirically identified wavelet temporal model was obtained using data collected at a 1 s sampling rate, some of the fast dynamics of the plant that gave rise to jumps were not recorded. The fast dynamics comes primarily from the spatially distributed riser geometry. Hence, we simulated the impulse response of the 2-PDE riser model [37], approximated the faster riser dynamics with the transient spatiotemporal NARMA-L1 [38] model, and used the result in a convolution to obtain a spatiotemporal model of the transients. We then put the empirical temporal NARX model and the fast dynamics spatiotemporal model in parallel. Finally, we combined the temporal GPC control and the spatiotemporal deadbeat control, each based on its respective model, into the closed loop dual-model self-tuning configuration shown in Figure 13. In this configuration, the dynamic inter-loop coupling is rather minimal due to significantly differing time scales of each elf-tuning loop, ideally requiring two-sampling-rates hardware/software setting, not available for this experiment.

The nonlinear 2-PDE model governing the evolution of the variables (voidage and solid velocity) in the riser can be represented [37] as:

$$\begin{array}{l}\frac{\partial \epsilon }{\partial t}=(1-\epsilon )\frac{\partial u_s}{\partial x}-u_s\frac{\partial \epsilon }{\partial x},\\ \frac{\partial u_s}{\partial t}=-u_s\frac{\partial u_s}{\partial x}+C_1{\epsilon }^{-6.7}-C_2{\epsilon }^{-5.7}{u}_s+C_3{\epsilon }^{-4.7}{u}_s^2+C_4{(1-\epsilon )}^{-0.54}-C_5,\end{array}$$

(55)

where $\epsilon$ is the voidage and $u_s$ is the solid velocity. The other parameters are defined in [37]. From simulations of this model, we could obtain a response $h(x,t)$ to an impulse actuation in solid velocity with area of 0.1.

Then the response to an arbitrary inlet solid velocity $u(t)$ can be calculated as:

$$y(x,t)={\displaystyle {\int }_{-\infty }^{t}h(x,\tau )}\cdot 10u(t-\tau )d\tau ,$$

(56)

where the scaling factor is necessary since the simulated input was not 1. The simulated impulse responses are shown in Figure 14.

Since the impulse response is uniformly zero after 0.6 s, Equation (56) can be limited to:

$$y(x,t)={\displaystyle {\int }_{t-0.6}^{t}h(x,\tau )}\cdot 10u(t-\tau )d\tau .$$

(57)

To obtain a low-order high fidelity finite-dimensional representation of the impulse response, a wavelet decomposition [39,40] was used to approximate $h(x,t)$. That is, the impulse response $h(x,t)$ was approximated as:

$$h(x,t)={\displaystyle \sum _{m=1}^{m_{\max }}{\displaystyle \sum _{n=1}^{n_{\max }}{\beta }_m(x)}}{c}_{m,n}{\alpha }_n(t),$$

(58)

where $\left\{{\beta }_m(x)\right\}$ and $\left\{{\alpha }_n(t)\right\}$ are wavelet basis functions.

Figure 15 is the resulting wavelet approximation of the impulse response. Here we chose Gaussian wavelet functions specifically. In this case, 23 spatial and 22 temporal wavelets were used. The coefficients $c_{m,n}$ were determined using a least-squares regression.

The following notation is used to divide the impulse response into separate parts for voidage $\epsilon$ and velocity $u_s$:

$$\begin{array}{l}{h}_{u_s}(x,t)={\displaystyle \sum _{m=1}^{m_{\max }}{\displaystyle \sum _{n=1}^{n_{\max }}{\beta }_m(x)}}{c}_{m,n}{\alpha }_n(t),\\ h_{\epsilon }(x,t)={\displaystyle \sum _{m=1}^{m_{\max }}{\displaystyle \sum _{n=1}^{n_{\max }}{\beta }_m(x)}}{d}_{m,n}{\alpha }_n(t).\end{array}$$

(59)

Then using the convolution:

$$\mathsf{\Delta }{u}_s(x,t)=h_{u_s}(x,\tau )\ast 10(u-\tau )=10{\displaystyle {\int }_0^{0.6}{\displaystyle \sum _{m=1}^{m_{\max }}{\displaystyle \sum _{n=1}^{n_{\max }}{\beta }_m(x)}}{c}_{m,n}{\alpha }_n(\tau )}u(t-\tau )d\tau .$$

(60)

Since the online measurements were only available at 1 s intervals, it was assumed that:

$$u(t-\tau )=(1-\tau )u(t)+\tau u(t-1), 0\le \tau \le 1,$$

(61)

to linearly interpolate between the measurements. Then:

$$\begin{array}{cc}\hfill \mathsf{\Delta }{u}_s(x,t)& =10{\displaystyle {\int }_0^{0.6}{\displaystyle \sum _{m=1}^{m_{\max }}{\displaystyle \sum _{n=1}^{n_{\max }}{\beta }_m}}{c}_{m,n}{\alpha }_n(\tau )\left[(1-\tau )u(t)+\tau u(t-1)\right]}d\tau \hfill \\ & =10{\displaystyle \sum _{m=1}^{m_{\max }}{\displaystyle \sum _{n=1}^{n_{\max }}{\beta }_m(x)}}{c}_{m,n}{\displaystyle {\int }_0^{0.6}{\alpha }_n(\tau )\left[(1-\tau )u(t)+\tau u(t-1)\right]}d\tau \hfill \\ & =10\left[{\displaystyle \sum _{m=1}^{m_{\max }}{\displaystyle \sum _{n=1}^{n_{\max }}{\beta }_m(x)}}{c}_{m,n}{\displaystyle {\int }_0^{0.6}(1-\tau ){\alpha }_n(\tau )}d\tau \right]u(t)+10\left[{\displaystyle \sum _{m=1}^{m_{\max }}{\displaystyle \sum _{n=1}^{n_{\max }}{\beta }_m(x)}}{c}_{m,n}{\displaystyle {\int }_0^{0.6}\tau {\alpha }_n(\tau )}d\tau \right]u(t-1).\hfill \end{array}$$

(62)

Denote:

$$\begin{array}{l}{\alpha }_{n,0}={\displaystyle {\int }_0^{0.6}(1-\tau )}{\alpha }_n(\tau )d\tau ,\\ {\alpha }_{n,1}={\displaystyle {\int }_0^{0.6}\tau }{\alpha }_n(\tau )d\tau ,\end{array}$$

(63)

and:

$$\begin{array}{l}{\gamma }_0(x)=10{\displaystyle \sum _{m=1}^{m_{\max }}{\displaystyle \sum _{n=1}^{n_{\max }}{\beta }_m(x)}{c}_{m,n}{\alpha }_{n,0}},\\ {\gamma }_1(x)=10{\displaystyle \sum _{m=1}^{m_{\max }}{\displaystyle \sum _{n=1}^{n_{\max }}{\beta }_m(x)}{c}_{m,n}{\alpha }_{n,1}}.\end{array}$$

(64)

Then Equation (62) simplifies to:

$$\mathsf{\Delta }{u}_s(x,t)={\gamma }_0(x)u(t)+{\gamma }_1(x)u(t-1).$$

(65)

Similarly, we have:

$$\mathsf{\Delta }\epsilon (x,t)={\eta }_0(x)u(t)+{\eta }_1(x)u(t-1),$$

(66)

where:

$$\begin{array}{l}{\eta }_0(x)=10{\displaystyle \sum _{m=1}^{m_{\max }}{\displaystyle \sum _{n=1}^{n_{\max }}{\beta }_m(x)}{d}_{m,n}{\alpha }_{n,0}},\\ {\eta }_1(x)=10{\displaystyle \sum _{m=1}^{m_{\max }}{\displaystyle \sum _{n=1}^{n_{\max }}{\beta }_m(x)}{d}_{m,n}{\alpha }_{n,1}}.\end{array}$$

(67)

Omitting for brevity several routine manipulations (available in [23]), the output $\mathrm{DP}47$ representing the pressure drop across the riser can be given as:

$$\begin{array}{ll}\hfill P(5)-P(0)=& \left[{\rho }_g\epsilon (0)u_g^2(0)+{\rho }_s(1-\epsilon (0))u_s^2(0)\right]\hfill \\ & -\left[{\rho }_g\epsilon (5)u_g^2(5)+{\rho }_s(1-\epsilon (5))u_s^2(5)\right]\hfill \\ & -{\displaystyle {\int }_0^{5}g\left[{\rho }_g\epsilon (y)+{\rho }_s(1-\epsilon (y))\right]}dy,\hfill \end{array}$$

(68)

where the riser length of 5 m is used, constant $g$ is the gravity acceleration, $u$ is the velocity, $\rho$ is the density, the subscripts $s$ and $g$ refer to solid and gas, respectively, and $u_g=\frac{U_g}{\epsilon (x)}$ where $U_g$ is the superficial gas velocity. Expanding Equation (68) yields:

$$\begin{array}{cc}\hfill P(5)-P(0)=& \left[{\rho }_g\frac{U_g^2}{\epsilon (0)}+{\rho }_s(1-\epsilon (0))u_s^2(0)\right]\hfill \\ & -\left[{\rho }_g\frac{U_g^2}{\epsilon (5)}+{\rho }_s(1-\epsilon (5))u_s^2(5)\right]\hfill \\ & -{\displaystyle {\int }_0^{5}g}\left[{\rho }_g\epsilon (y)+{\rho }_s(1-\epsilon (y))\right]dy\hfill \\ & ={\rho }_g\frac{U_g^2}{{\epsilon }_{ss}(0)+\mathsf{\Delta }\epsilon (0)}+{\rho }_s(1-{\epsilon }_{ss}(0)-\mathsf{\Delta }\epsilon (0)){(u_s(0)+\mathsf{\Delta }{u}_s(0))}^2\hfill \\ & -{\rho }_g\frac{U_g^2}{{\epsilon }_{ss}(5)+\mathsf{\Delta }\epsilon (5)}-{\rho }_s(1-{\epsilon }_{ss}(5)-\mathsf{\Delta }\epsilon (5)){(u_s(5)+\mathsf{\Delta }{u}_s(5))}^2\hfill \\ & -g{\rho }_g{\displaystyle {\int }_0^5({\epsilon }_{ss}(y)+\mathsf{\Delta }\epsilon (y))}dy-g{\rho }_s{\displaystyle {\int }_0^5(1-{\epsilon }_{ss}(y)-\mathsf{\Delta }\epsilon (y))}dy,\hfill \end{array}$$

(69)

where the subscript $ss$ designates the steady state. Now, substituting the wavelet model gives:

$$\begin{array}{cc}\hfill P(5)-P(0)=& {\rho }_g\frac{U_g^2}{{\epsilon }_{ss}(0)+{\eta }_0(0)u(t)+{\eta }_1(0)u(t-1)}\hfill \\ & +{\rho }_s\left[1-{\epsilon }_{ss}(0)-{\eta }_0(0)u(t)-{\eta }_1(0)u(t-1)\right]{\left[u_s(0)+{\gamma }_0(0)u(t)+{\gamma }_1(0)u(t-1)\right]}^2\hfill \\ & -{\rho }_g\frac{U_g^2}{{\epsilon }_{ss}(5)+{\eta }_0(5)u(t)+{\eta }_1(5)u(t-1)}\hfill \\ & -{\rho }_s\left[1-{\epsilon }_{ss}(5)-{\eta }_0(5)u(t)-{\eta }_1(5)u(t-1)\right]{\left[u_s(5)+{\gamma }_0(5)u(t)+{\gamma }_1(5)u(t-1)\right]}^2\hfill \\ & -g{\rho }_g{\displaystyle {\int }_0^5\left({\epsilon }_{ss}(y)+{\eta }_0(y)u(t)+{\eta }_1(y)u(t-1)\right)}dy\hfill \\ & -g{\rho }_s{\displaystyle {\int }_0^5\left(1-{\epsilon }_{ss}(y)-{\eta }_0(y)u(t)-{\eta }_1(y)u(t-1)\right)}dy.\hfill \end{array}$$

(70)

Our goal was to use the model given by Equation (70) to account for the spatiotemporal behavior of the CL system to the extent allowed by the available sampling rate, and also to develop the control setting to be ready to employ much higher sampling rates for performance improvement, once they become available on the test rig through the processor upgrades to the GPUs and FPGAs. It is also of interest to calculate the steady-state pressure drop:

$$\begin{array}{cc}\hfill \mathsf{\Delta }{P}_0=& {\rho }_g\frac{U_g^2}{{\epsilon }_{ss}(0)+{\eta }_0(0)u(t)+{\eta }_1(0)u(t)}\hfill \\ & +{\rho }_s\left[1-{\epsilon }_{ss}(0)-{\eta }_0(0)u(t)-{\eta }_1(0)u(t)\right]{\left[u_s(0)+{\gamma }_0(0)u(t)+{\gamma }_1(0)u(t)\right]}^2\hfill \\ & -{\rho }_g\frac{U_g^2}{{\epsilon }_{ss}(5)+{\eta }_0(5)u(t)+{\eta }_1(5)u(t)}\hfill \\ & -{\rho }_s\left[1-{\epsilon }_{ss}(5)-{\eta }_0(5)u(t)-{\eta }_1(5)u(t)\right]{\left[u_s(5)+{\gamma }_0(5)u(t)+{\gamma }_1(5)u(t)\right]}^2\hfill \\ & -g{\rho }_g{\displaystyle {\int }_0^5\left({\epsilon }_{ss}(y)+{\eta }_0(y)u(t)+{\eta }_1(y)u(t)\right)}dy\hfill \\ & -g{\rho }_s{\displaystyle {\int }_0^5\left(1-{\epsilon }_{ss}(y)-{\eta }_0(y)u(t)-{\eta }_1(y)u(t)\right)}dy.\hfill \end{array}$$

(71)

This is the pressure drop predicted by this model for constant input $u\left(t\right)$, as opposed to the linear interpolation described above. We can then use this model to approximate the transient difference, and the NARX wavelet model to approximate the steady state. The difference between the transient pressure drop $\mathsf{\Delta }P(t)$ and the eventual steady pressure drop $\mathsf{\Delta }{P}_0(t)$ is then equal to:

$$\mathsf{\Delta }P-\mathsf{\Delta }{P}_0\hspace{0.17em}=\hspace{0.17em}P(5)-P(0)-\mathsf{\Delta }{P}_0.$$

(72)

Linearizing Equation (72) about $u(t)=u(t-1)$ gives:

$$\mathsf{\Delta }P-\mathsf{\Delta }{P}_0\approx f\left(u(t-1)\right)\left(u(t)-u(t-1)\right),$$

(73)

where:

$$\begin{array}{cc}\hfill f(y)=& k_1\frac{(k_2+k_{3}y)}{{(k_2+k_{4}y)}^4}+k_5{(k_6+k_{7}y)}^2\hfill \\ & +k_8(1-k_2-k_{4}y)(k_6+k_{9}y)+k_{10}\frac{(k_{11}+k_{12}y)}{{(k_{11}+k_{13}y)}^4}\hfill \\ & +k_{14}{(k_{15}+k_{16}y)}^2+k_{17}(1-k_{11}-k_{13}y)(k_{15}+k_{18}y)+k_{19},\hfill \end{array}$$

(74)

and:

$$\begin{array}{l}{k}_1={\rho }_g{U}_g{\eta }_1(0), k_2={\epsilon }_{ss}(0), k_3={\eta }_0(0)+2{\eta }_1(0),\\ k_4={\eta }_0(0)+{\eta }_1(0), k_5={\rho }_s{\eta }_1(0), k_6=u_{s,ss}(0),\\ k_7={\gamma }_0(0)+{\gamma }_1(0), k_8=-{\rho }_s{\gamma }_1(0), k_9={\gamma }_0(0)+2{\gamma }_1(0),\\ k_{10}=-{\rho }_g{U}_g{\eta }_1(5), k_{11}={\epsilon }_{ss}(5), k_{12}={\eta }_0(5)+2{\eta }_1(5),\\ k_{13}={\eta }_0(5)+{\eta }_1(5), k_{14}=-{\rho }_s{\eta }_1(5), k_{15}=u_{s,ss}(5),\\ k_{16}={\gamma }_0(5)+{\gamma }_1(5), k_{17}={\rho }_s{\gamma }_1(5), k_{18}={\gamma }_0(5)+2{\gamma }_1(5),\\ k_{19}=g({\rho }_g-{\rho }_s){\displaystyle {\int }_0^5{\eta }_1(y)}\hspace{0.17em}dy.\end{array}$$

(75)

The input to the computational model was in terms of the velocity boundary condition, so that $u(t)=\mathsf{\Delta }{u}_s(0,t)$. This can be connected to the inputs $S_1$ and $S_2$ via a quadratic model fitted to the test data where:

$$\begin{array}{ll}\hfill u(t)\approx & \frac{1}{{\epsilon }_0}\left(2a_1{S}_1(t-1)+a_3{S}_2(t)+a_4\right)\left(S_1(t)-S_1(t-1)\right)\hfill \\ & +\frac{1}{{\epsilon }_0}\left(a_1{S}_1^2(t-1)+a_2{S}_2^2(t)+a_3{S}_1(t-1)S_2(t)+a_4{S}_1(t-1)+a_5{S}_2(t)+a_6\right)-u_{s,ss}(0)\hfill \\ & =\frac{1}{{\epsilon }_0}\left(a_1{S}_1(t-1)+a_3{S}_2(t)+a_4\right)S_1(t)\hfill \\ & +\frac{1}{{\epsilon }_0}\left(a_2{S}_2^2(t)+a_5{S}_2(t)+a_6\right)-u_{s,ss}(0).\hfill \end{array}$$

(76)

Then:

$$\begin{array}{ll}\hfill \mathsf{\Delta }P-\mathsf{\Delta }{P}_0& \approx \frac{1}{{\epsilon }_0}f\left(u(t-1)\right)\left(2a_1{S}_1(t-1)+a_3{S}_2(t)+a_4\right)S_1(t)\hfill \\ & +f(u(t-1))\left[\frac{1}{{\epsilon }_0}\left(a_2{S}_2^2(t)+a_5{S}_2(t)+a_6\right)-u_{s,ss}(0)-u(t-1)\right]\hfill \\ & =g_{\mathsf{\Delta }P}\left(S_1(t-1),S_2(t-1),S_2(t)\right)S_1(t)+f_{\mathsf{\Delta }P}\left(S_1(t-1),S_2(t-1),S_2(t)\right).\hfill \end{array}$$

(77)

The NARX wavelet MRA model takes the form as in Equation (51):

$$y_w(t)=f\left(y(t-1),\cdots ,y(t-n_y),S_{1\_w}(t-1),\cdots ,S_{1\_w}(t-n_u)\right),$$

(78)

where $S_{1,w}(t)$ is the control command calculated by wavelet adaptive GPC control. Then, the fast transient behavior model of Equation (77) can be combined with Equation (78) to obtain a spatiotemporal multiscale dynamic network representation of the entire CL process:

$$y(t)\approx y_w(t)-f_{\mathsf{\Delta }P}-g_{\mathsf{\Delta }P}{S}_1(t).$$

(79)

The sign change is necessary because the pressure drop across the riser is negative in the model above, i.e., $P(5)-P(0)<0$. Then, the deadbeat predictive controller taking account of fast dynamics is:

$$S_{1,\mathrm{fast}}(t)=\frac{y_r(t)-y_w(t)+f_{\mathsf{\Delta }P}}{g_{\mathsf{\Delta }P}},$$

(80)

where $y_r(t)$ is the reference signal. Hence, the final spatiotemporal wavelet controller $S_1(t)$ implemented on the real CL process is given by:

$$S_1(t)=S_{1,w}(t)+S_{1,\mathrm{fast}}(t).$$

(81)

This controller was also implemented in the single loop cold flow CL test rig. $S_1$ was taken as the single input and $\mathrm{DP}47$ as the output, while $S_2$ was mostly steady, but with jumps. The reference signal was set to 16 initially and then reduced to 13 around time 17:01. The tracking response of system output and the corresponding control efforts are shown in Figure 16 and Figure 17, respectively. The controller is seen to stabilize the system quite well under difficult operating conditions. The pressure drop DP47 over the riser related to the fluidizing air flows was discussed based on the closed loop topology augmentation with the spatiotemporal model-based control.

5. Discussion

The control objective of developing the tracking controller for the CL test rig at ALSTOM Power was addressed at the start of the project through both first-principles model development and empirical system identification from the input/output experimental data record. The first approach resulted in the analytical development and the numerical simulation of the 1D, 2D, and 3D partial differential equation (PDE) networks—systems of coupled PDEs, each describing the testbed subsystem. On the basis of these models, approximately tuned to the process through the experimental data, first, the linear finite dimensional models were developed, after which robust controllers based on the $H_{\infty }$ approach were designed. The latter were implemented on the test rig; however, no satisfactory performance was obtained. The empirical approach initially involved a polynomial NARX model test fit with the use of model predictive control (MPC) with constraints attained through the QP (quadratic programming) based control signal calculation. At the same time the hypothesis was posed that the process is highly multiscale and that a multiscale controller design should be attempted. The experimental data was then subjected to multi-resolution decomposition and was indeed found to be highly multiscale [23]. After some debate, it was suggested to refocus the effort from model-based robust controller improvement and traditional constrained MPC to empirical multiresolution controller development. The initial step was to fit a multiresolution nonlinear, but linear in parameters, temporal model to the input/output data, with the subsequent steps involving the self-tuning GPC-type controller development, as shown in Figure 3. After a number of prolonged experiments, the resulting controller was outfitted with rate constraints and tuned to work on the test rig, showing performance noticeably superior to that of the other techniques employed. Further on, since the underlying process dynamics was exhibited by the intrinsically spatially distributed processes, a spatiotemporal multiresolution control component was developed for the riser dynamics and combined with the temporal topology. The latter part was implemented and tested, producing reasonable overall performance, but it could not be fully utilized due to the computational real-time controller limitations. Subsequently, the results presented in this paper were acknowledged by ALSTOM as making a real breakthrough in the project. The findings and their implications in the broad context imply that if a robust linear controller does not work well on a system, the system is likely rather complex, and the nonlinear multiresolution modeling in combination with linear-type control structures under constraints could offer effective configurations for high performance system control.

The future theoretical effort in the temporal implementation in Section 3 should involve developing the proof of the simultaneous error convergence of the coupled identifier/controller system under control rate constraints—a challenging analytical task. Future research should also address spatiotemporal controller development and implementation in more generality than that presented in Section 4. In the present work, the latter controller had limited efficacy due to the low sampling/computation rate of the test rig control processor, but the approach presented stands ready for further theoretical advancement and applications involving more advanced hardware.

Due to the combination of the coarser temporal and the finer spatiotemporal control in the overall control configuration, the results presented have a broad appeal in other applications involving complex multiscale spatiotemporal dynamics, such as, for example, robotic electrosurgery [41]. In this application, the motion of the cutting probe is strongly coupled to the spatiotemporal electrosurgical impact of the latter on the target tissues, and the near-field probe-tissue interaction process is best described by a complex time-fractional PDE [42]. These features make the technique proposed a good match for this application.

6. Conclusions

In this paper, the following results have been presented:

• closing the gap between the actual system output data and that predicted by the first-principles model of the complex chemical looping process through the empirically identified wavelet multiresolution analysis model initially trained on-line to estimate the nonlinear dynamic characteristics;
• embedding the multiresolution model into the generalized predictive control structure to obtain self-tuning control strategy for stable tracking of a chemical looping process with complex process dynamics under actuator rate constraints;
• showing boundedness of the adaptation gains for identification and control laws using the Lyapunov function theorems, and providing a guidance for attainment of asymptotic stability of the closed-loop system through the choice of these adaptation gains;
• using a spatiotemporal wavelet decomposition of the impulse response of the chemical looping process riser for designing a deadbeat predictive controller for further enhancement of the closed loop performance to account for fast system dynamics.
• experimentally confirming the effectiveness of the proposed controller design methods though their implementation on the novel single loop chemical looping cold flow testbed with complex dynamics,

Limitations of the present study lie in the insufficient spatiotemporal modeling and controller design for the riser and the inability to fully utilize the designed spatiotemporal controller for it because of the insufficient real-time performance of the control processor and the data acquisition system. Future efforts should focus on the advancement of the spatiotemporal part of the design, the overall controller robustness evaluation and enhancement, the development of the rigorous convergence proofs for the coupled identification/control laws under control rate constraints for the temporal and the entire spatiotemporal control topologies, and the implementation of the controller using advanced processors. The proposed techniques is planned to be applied by the authors to other areas, such as robotic electrosurgery.