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 varyingH∞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

_ϕj,k(x): orthonormal basis for_Vj;_ψj,k(x): orthonormal basis for_Wj;y(t): system output;u(t): system input;w(t): system noise;e(t): model output estimation error;_ny: maximum lag in the output;_nu: maximum lag in the input;θ: weighting parameter vector trained on-line;_gi: multivariable scaling or wavelet basis function of past inputs and outputs;_γθ: adaptation gain for the control input vector;U: control input vector;Δ_utarget: unconstraint control signal calculated by the predictive control law;μ: design parameter;ε:voidage;_us: solid velocity;_Ug: superficial gas velocity;_S1,w(t): control command calculated by wavelet adaptive GPC control;_yr(t): reference signal;_{∇θ J1}(n): gradient of loss function_J1with respect toθ;_{∇U J2}(n): gradient of the loss function_J2with respect toU;_∇θe(n): sensitivity derivative at timen.

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) includeS1,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 isDP47, 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_Vj∈_L2(R),j∈Zsatisfying:

⋯_V−1⊂_V0⊂_V1⋯,

with the following properties:

_{∪j∈Z Vj} is dense in _L2(R); _{∩j∈Z Vj}={0},

f(x)∈_Vj⇔f(2x)∈_Vj+1

f(x)∈_Vj⇔f(x−^2−jk)∈_Vj, k∈Z,

_Vj=span{_ϕj,k,k∈Z},

whereZis the set of all integers,_ϕj,k(x)=^2j/2ϕ(^2jx−k)is an orthonormal basis for_Vjand_L2(R)is the space of square integrable functions of scalar real variable.

If we define_Wjto be the orthogonal complement of_Vjin_Vj+1, then:

_Vj+1=_Vj⊕_Wj, _Vj⊥_Wj,

_Wj=span{_ψj,k,k∈Z},

where_ψj,k(x)=^2j/2ψ(^2jx−k)is an orthonormal basis for_Wj. It follows from Equations (1) and (6) that, any_Vjcan be written for anyl<jas:

_Vj+1=_Vl⊕_Wl⊕_Wl+1⊕_Wl+2⊕⋯⊕_Wj,

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

_L2(R)=_{⊕j∈Z Wj}

The functions_ϕj,kand_ψj,kwill be referred to as scaling and wavelet functions respectively. According to Equations (8) and (9), anyf(x)∈_L2(R)can be represented as:

f(x)=∑_n(f,_ϕJ,n)_ϕJ,n+∑_j≥J,n(f,_ψj,n)_ψj,n

The approximation starts from some lower resolution levelJand can be truncated at certain higher resolution levelNwhen:

∥f(x)−[∑n(f,_ϕJ,n)_ϕJ,n(x)+∑j=JN∑n(f,_ψj,n)_ψj,n(x)]∥<ε,

for any predefined small errorε>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),⋯,y(t−_ny),u(t−1),⋯,u(t−_nu))+w(t),

wherey(t),u(t),w(t)are the system output, input, and noise, andtis discrete time, respectively,_nyand_nuare the maximum lags in the output and input,w(t)is assumed to be a zero mean, independent, and bounded noise variable, andf(⋅)is some nonlinear function. Unless some prior knowledge of the system dynamics is available, most methods use nonparametric regression to estimate the nonlinear functionffrom the data. In our case,fis implemented as a linear expansion in terms of the scaling and wavelet functions of regressors_gisuch that

f=_{∑i=1m θi gi}

minimizes a pre-specified approximation adequacy criterion, whereθ={_θi}is a parameter vector trained on-line,_gi∈{_ϕj,k,_ψj,k}is a multivariable scaling or wavelet basis function of past inputs and outputs, andmis 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, letDP47be the actual system outputyandS1be the control inputu, whileS2is set to a constant value. Lety^denote the approximated system output. Then, the identified wavelet MRA based model is defined as follows:

y^(t)=f(y(t−1),⋯,y(t−_ny),u(t−1),⋯,u(t−_nu)),

where f is defined in Equation (13). Then, the error between the real plant outputyand the estimated outputy^is defined as:

e(n)=y(n)−y^(n).

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

_J1(n)=12^e2(n),

wherenindicates discrete time. To make_J1small, we employ a parameter adaptation law in the form of a gradient descent (GD) algorithm, which adjusts the weighting gainsθto keep the gradient of_J1negative, that is:

θ(n+1)=θ(n)−_{γθ ∇θ J1}(n)=θ(n)−_γθe(n)_∇θe(n),

where_γθis the adaptation gain,_{∇θ J1}(n)is the gradient of_J1with respect toθat discrete timen, and_∇θe(n)is the so-called sensitivity derivative at timenindicating how the error is influenced by the weighting parametersθ. From Equations (13)–(15), the sensitivity derivative_∇θecan be derived as follows:

_∇θe=−_∇θy^=−_∇θf=−g⇒θ(n+1)=θ(n)+_γθe(n)g(n).

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

_J2=12{_∑k=N1N2 ^{(_ym(n+k)−y^(n+k))2}+_{∑k=1Nu ρk}Δu^(n+k−1)2},

where_N1and_N2are the minimum and the maximum output prediction horizons, respectively,_Nuis the control horizon,Δis the difference operator,Δu(n)=u(n)−u(n−1), and_ρkis thek-th control weighting factor. Assuming_N1=1,_N2=L=_Np, and identical control weighing factor_ρk=ρ,Equation (19) can be rewritten in the vector form as:

_J2=12{^{∥_Ym(n+1)−Y^(n+1)∥2}+ρ^{∥ΔU(n)∥2}},

where:

_Ym(n+1)=^{[_ym(n+1),_ym(n+2),⋯,_ym(n+L)]T},

Y^(n+1)=^{[y^(n+1),y^(n+2),⋯,y^(n+L)]T},

U(n)=^{[u(n),u(n+1),⋯,u(n+_Nu−1)]T},

ΔU(n)=^{[Δu(n),Δu(n+1),⋯,Δu(n+_Nu−1)]T},

and∥⋅∥is the_L2norm of then-dimensional real vectors.

The loss function_J2is now minimized to drive the system outputyto the reference signal_ymgiven 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_J2, the GD method is implemented again to recursively calculate the_Nu-dimensional control increment sequenceΔUas follows:

ΔU(n)=−_{γu ∇U J2}(n),

where_γuis the adaptation gain for the control input vectorU. Noting that for any vector y(x),_∇x ^∥y∥2=2(_∇xy)y, from Equations (20) and (21), the gradient of the loss function_J2with respect toUcan be obtained as:

_{∇U J2}(n)=12_∇U{^{∥_Ym(n+1)−Y^(n+1)∥2}+ρ^{∥ΔU(n)∥2}}=_∇U{_Ym(n+1)−Y^(n+1)}{_Ym(n+1)−Y^(n+1)}+ρ_∇U{ΔU(n)}{ΔU(n)}.

The first part of the expression above is evaluated as follows. First, we note that since_ym(n)are the reference signals, which are preset constants, we have∂_ym(n+1)/∂u(n)=0. Then, sincey^(n+k)=^θTu,y^depends only on the past u, we have:

∂y(n+k)∂u(n+l)={∂y^(n+k)∂u(n+l), when k>l,0, when k≤l.

This yields:

_∇U{_Ym(n+1)−Y^(n+1)}=

_{[∂∂u(n){ym(n+1)−y^(n+1)}∂∂u(n){ym(n+2)−y^(n+2)}⋯∂∂u(n){ym(n+L)−y^(n+L)}∂∂u(n+1){ym(n+1)−y^(n+1)}∂∂u(n+1){ym(n+2)−y^(n+2)}⋯∂∂u(n+1){ym(n+L)−y^(n+L)}⋮⋮⋱⋮∂∂u(n+Nu−1){ym(n+1)−y^(n+1)}∂∂u(n+Nu−1){ym(n+2)−y^(n+2)}⋯∂∂u(n+Nu−1){ym(n+L)−y^(n+L)}]Nu×Np}

=_{[−∂y^(n+1)∂u(n)−∂y^(n+2)∂u(n)⋯−∂y^(n+Nu)∂u(n)⋯−∂y^(n+L)∂u(n)0−∂y^(n+2)∂u(n+1)⋯−∂y^(n+Nu)∂u(n+1)⋯−∂y^(n+L)∂u(n+1)⋮0⋯0−∂y^(n+Nu)∂u(n+Nu−1)⋯−∂y^(n+L)∂u(n+Nu−1)]Nu×Np}

The second part of_{∇U J2}(n)is evaluated by taking into account the relationΔu(n)=u(n)−u(n−1),so that∂Δu(n)/∂u(n)=1and∂Δu(n)/∂u(n−1)=−1. The latter yields:

_∇U{ΔU(n)}=_{[∂Δu(n)∂u(n)∂Δu(n+1)∂u(n)∂Δu(n+2)∂u(n)⋯∂Δu(n+Nu−1)∂u(n)∂Δu(n)∂u(n+1)∂Δu(n+1)∂u(n+1)∂Δu(n+2)∂u(n+1)⋯∂Δu(n+Nu−1)∂u(n+1)⋮⋮⋮⋱⋮∂Δu(n)∂u(n+Nu−1)⋯∂Δu(n+Nu−3)∂u(n+Nu−1)∂Δu(n+Nu−2)∂u(n+Nu−1)∂Δu(n+Nu−1)∂u(n+Nu−1)]Nu×Nu}

=_{[1−10⋯001−1⋯0⋮⋮⋮⋮⋮0⋯001]Nu×Nu}.

Combining the above expressions into:

_{∇U J2}(n)=−G(_Ym(n+1)−Y^(n+1))+ρHΔU(n),

and substituting into Equation (21) as:

ΔU(n)=−_{γu ∇U J2}(n)=_γuG(_Ym(n+1)−Y^(n+1))−_γuρHΔU(n),

where:

G=_{[∂y^(n+1)∂u(n)∂y^(n+2)∂u(n)⋯∂y^(n+Nu)∂u(n)⋯∂y^(n+L)∂u(n)0∂y^(n+2)∂u(n+1)⋯∂y^(n+Nu)∂u(n+1)⋯∂y^(n+L)∂u(n+1)⋮0⋯0∂y^(n+Nu)∂u(n+Nu−1)⋯∂y^(n+L)∂u(n+Nu−1)]Nu×Np}

and:

H=_{[1−10⋯001−1⋯0⋮⋮⋮⋮⋮0⋯001]Nu×Nu}

yields the control law of the form:

ΔU(n)=^{(I+_γuρH)−1} _γuG(_Ym(n+1)−Y^(n+1)),

whereIis the_Nu×_Nuidentity 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:

_V1(n)=12^e2(n),

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

Δ_V1(n)=_V1(n+1)−_V1(n)=12(^e2(n+1)−^e2(n)).

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

Δe(n)=e(n+1)−e(n)=[∂e(n)∂θ(n)]Δθ(n)

whereΔθ(n)=_{{Δθi(n)}i=1m}represents a change in the arbitrary component of the weighting gain vectorθ. From Equation (18),Δ_θi(n)can be obtained by:

Δθ(n)=_γθe(n)g(n),

[∂e(n)∂θ(n)]=−[∂y^(n)∂θ(n)]=−^gT(n),

whereg(n)=_{gi(n)}i=1m.

Theorem 1.

Let_γθbe the adaptation gain for the weights of the wavelet MRA identified model and_gmaxbe defined as_gmax:=_maxn∥g(n)∥, wheregis the wavelet MRA basis function,nis the discrete time index, and∥⋅∥is the_L2norm of a real vector. Then convergence is guaranteed if_γθis chosen as:

0<_γθ<2_gmax2

Proof.

From Equations (28)–(31),Δ_V1(n)can be represented as:

Δ_V1(n)=Δe(n)[e(n)+12Δe(n)]=[∂e(n)∂θ]_γθe(n)g(n)×{e(n)+12[∂e(n)∂θ]_γθe(n)g(n)}

=−_γθ ^e2(n)^∥g(n)∥2+12_γθ2 ^e2(n)^∥g(n)∥4=−λ^e2(n),

where:

λ=12_γθ ^∥g(n)∥2(2−_γθ ^∥g(n)∥2)≥12_γθ ^∥g(n)∥2(2−_{γθ gmax2})>0.

From Equation (32) we obtain_V1(n)≥0andΔ_V1(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:

_V2(n)=12^{∥E(n+1)∥2},

whereE(n+1)=_Ym(n+1)−Y^(n+1). Then the change of the Lyapunov function is obtained as:

Δ_V2(n)=_V2(n+1)−_V2(n)=12(^{∥E(n+2)∥2}−^{∥E(n+1)∥2}).

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

ΔE(n+1)=E(n+2)−E(n+1)=[∂E(n+1)∂U(n)]ΔU(n),

whereΔU(n)is defined in Equation (26) and∂E(n+1)∂U(n)=−^GT. Then Equation (37) can be expressed as:

ΔE(n+1)=−^{GT (I+_γuρH)−1} _γuGE(n+1).

Theorem 2.

Let_γube the adaptation gain for the GPC control input sequence. Assume a control weighting factorρ>0. Then the stable tracking convergence of the wavelet MRA based GPC control system is guaranteed if:

0<_γu<2_λmax(G^GT),

where_λmax(⋅)is the maximum eigenvalue of the matrix.

Proof.

From Equations (36)–(38),Δ_V2(n)can be represented as:

Δ_V2(n)=12[^{(E(n+1)+ΔE(n+1))T}(E(n+1)+ΔE(n+1))−E^(n+1)TE(n+1)]=Δ^ET(n+1)[E(n+1)+12ΔE(n+1)]=−^(GE)T _γu ^{((I+_γuρH)−1)T}[I−12G^{GT (I+_γuρH)−1} _γu]GE=−^(GE)T _{R1 R2}GE,

where:

_R1=_γu ^{((I+_γuρH)−1)T},

_R2=I−12_γuG^{GT (I+_γuρH)−1}.

☐

If_R1and_R2are both positive definite matrices, it follows thatΔ_V2(n)<0. Together with_V2(n)>0, the stable tracking of the reference signals is guaranteed.

From Equation (25) it can be shown that the eigenvalues ofHare_λH=_{{1,⋯,1}Nu×1}. Then the eigenvalues of_R1can be derived as:

_λR1=_{{γu ^{(1+γuρ)−1},⋯,γu ^{(1+γuρ)−1}}Nu×1}

Hence, all eigenvalues of_R1are positive if_γu>0. It follows that_R1>0.

If0<_γu<2_λmax(G^GT), then:

I−12_γuG^GT>0.

From Equation (25) we have:

_γuρH>0.

Then from Equations (44) and (45)

I−12_γuG^GT+_γuρH>0.

Similarly to the way it was done for Equation (43), we can prove thatI+_γuρH>0. Then Equation (46) can be rewritten as:

(I+_γuρH)(I−12_γuG^{GT (I+_γuρH)−1})>0.

SinceI+_γuρH>0,(I−12_γuG^{GT (I+_γuρH)−1})>0follows. Now we have_V2(n)≥0andΔ_V2(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:

Δ_umin≤u(n+k)−u(n+k−1)≤Δ_umax,

_umin≤u(n+k)≤_umax,

where0≤k≤_Nu−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 inputuis subject to an input magnitude constraint saturation:

^u*(n)=sat[u(n)]={_umin if u(n)<_uminumax if u(n)>_umaxu(n) otherwise.

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 outputywas selected to be the riser pressure drop DP47 (inch H₂O). Fluidizing air flowS1(standard cubic feet per hour, scfh), was used as the single control input u, while the other air flowS2(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:

y^(t)=f(y(t−1),⋯,y(t−_ny),u(t−1),⋯,u(t−_nu))=∑i=1m_{θi gi},

wherefis the unknown nonlinear mapping to be identified,u(t)andy(t)are the sampled input and output sequences,_nyand_nuare the maximum lags in the output and the input to be determined, respectively;θ={_θi}is the parameter vector trained on-line,_gi∈{_ϕj,k,_ψj,k}is a multivariable scaling or wavelet basis function of past inputs and outputs, andmis 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_S1was generated in the form of a pseudo random binary signal (PRBS) changed about a nominal value, and the pressure drop across the riserDP47 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 betweenS1andDP47based on the experimental data. The regressor set was specified as:

y(t−1),y(t−2),u(t−1),u(t−2),⋯,u(t−8).

Hence_ny=2, _nu=8 . For MRA model we chose radial Marr scaling and wavelet functions [12]:

ϕ(x)=exp(−0.5^∥x∥2), ψ(x)=(dim(x)−^∥x∥2)exp(−0.5^∥x∥2).

The initial coarse layer indexJwas chosen to be 3, with the number of basis functions doubling when resolution was incremented by 1 starting with 10. The final resolution adopted wasK=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 inputuwas subjected to rate constraints of the form:

Δu=Δ_utarget×exp(1−μ∥Δ_utarget∥),

whereΔ_utargetis the unconstraint control signal calculated by the predictive control law andμ>0is 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 systemyasymptotically tracks the reference signal_ym. The cost function to be minimized is defined in (19). The design parameters for GPC configuration were chosen as_N1=1, _N2=10, _Nu=8, ρ=1. The adaptation gains derived from Theorems 1 and 2, were chosen as_γθ=0.01, _γu=0.1. The system was initially stable around level of_y0=13inch H₂O. Two setpoint step change experiments were performed consecutively. After 5 min, the setpoint was first increased to16inch H₂O and stayed at the latter value for about 7 min. Then it went back to the original level of13inch H₂O. The air flowS2 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 inDP47 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_ym(t)=13+2sin(2π×0.01×t), whileS2 (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:

∂ε∂t=(1−ε)∂_us∂x−_us∂ε∂x,∂_us∂t=−_us∂_us∂x+_C1 ^ε−6.7−_C2 ^ε−5.7 _us+_C3 ^ε−4.7 _us2+_C4 ^{(1−ε)−0.54}−_C5,

whereεis the voidage and_us is the solid velocity. The other parameters are defined in [37]. From simulations of this model, we could obtain a responseh(x,t)to an impulse actuation in solid velocity with area of 0.1.

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

y(x,t)=_∫−∞th(x,τ)⋅10u(t−τ)dτ,

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)=_∫t−0.6th(x,τ)⋅10u(t−τ)dτ.

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

h(x,t)=∑m=1_mmax∑n=1_nmaxβm(x)_{cm,n αn}(t),

where{_βm(x)}and{_αn(t)}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_cm,nwere determined using a least-squares regression.

The following notation is used to divide the impulse response into separate parts for voidageεand velocity_us:

_hus(x,t)=∑m=1_mmax∑n=1_nmaxβm(x)_{cm,n αn}(t),_hε(x,t)=∑m=1_mmax∑n=1_nmaxβm(x)_{dm,n αn}(t).

Then using the convolution:

Δ_us(x,t)=_hus(x,τ)∗10(u−τ)=10_∫00.6∑m=1_mmax∑n=1_nmaxβm(x)_{cm,n αn}(τ)u(t−τ)dτ.

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

u(t−τ)=(1−τ)u(t)+τu(t−1), 0≤τ≤1,

to linearly interpolate between the measurements. Then:

Δ_us(x,t)=10_∫00.6∑m=1_mmax∑n=1_{nmaxβmcm,n αn}(τ)[(1−τ)u(t)+τu(t−1)]dτ=10∑m=1_mmax∑n=1_nmaxβm(x)_{cm,n∫00.6αn}(τ)[(1−τ)u(t)+τu(t−1)]dτ=10[∑m=1_mmax∑n=1_nmaxβm(x)_cm,n∫00.6(1−τ)_αn(τ)dτ]u(t)+10[∑m=1_mmax∑n=1_nmaxβm(x)_cm,n∫00.6τ_αn(τ)dτ]u(t−1).

Denote:

_αn,0=_∫00.6(1−τ)_αn(τ)dτ,_αn,1=_∫00.6τ_αn(τ)dτ,

and:

_γ0(x)=10∑m=1_mmax∑n=1_nmaxβm(x)_{cm,n αn,0},_γ1(x)=10∑m=1_mmax∑n=1_nmaxβm(x)_{cm,n αn,1}.

Then Equation (62) simplifies to:

Δ_us(x,t)=_γ0(x)u(t)+_γ1(x)u(t−1).

Similarly, we have:

Δε(x,t)=_η0(x)u(t)+_η1(x)u(t−1),

where:

_η0(x)=10∑m=1_mmax∑n=1_nmaxβm(x)_{dm,n αn,0},_η1(x)=10∑m=1_mmax∑n=1_nmaxβm(x)_{dm,n αn,1}.

Omitting for brevity several routine manipulations (available in [23]), the outputDP47representing the pressure drop across the riser can be given as:

P(5)−P(0)=[_ρgε(0)_ug2(0)+_ρs(1−ε(0))_us2(0)]−[_ρgε(5)_ug2(5)+_ρs(1−ε(5))_us2(5)]−_∫05g[_ρgε(y)+_ρs(1−ε(y))]dy,

where the riser length of 5 m is used, constantgis the gravity acceleration,uis the velocity,ρis the density, the subscriptssandgrefer to solid and gas, respectively, and_ug=_Ugε(x)where_Ugis the superficial gas velocity. Expanding Equation (68) yields:

P(5)−P(0)=[_ρgUg2ε(0)+_ρs(1−ε(0))_us2(0)]−[_ρgUg2ε(5)+_ρs(1−ε(5))_us2(5)]−_∫05g[_ρgε(y)+_ρs(1−ε(y))]dy=_ρgUg2εss(0)+Δε(0)+_ρs(1−_εss(0)−Δε(0))^{(_us(0)+Δ_us(0))2}−_ρgUg2εss(5)+Δε(5)−_ρs(1−_εss(5)−Δε(5))^{(_us(5)+Δ_us(5))2}−g_ρg∫05(_εss(y)+Δε(y))dy−g_ρs∫05(1−_εss(y)−Δε(y))dy,

where the subscriptssdesignates the steady state. Now, substituting the wavelet model gives:

P(5)−P(0)=_ρgUg2εss(0)+_η0(0)u(t)+_η1(0)u(t−1)+_ρs[1−_εss(0)−_η0(0)u(t)−_η1(0)u(t−1)]^{[_us(0)+_γ0(0)u(t)+_γ1(0)u(t−1)]2}−_ρgUg2εss(5)+_η0(5)u(t)+_η1(5)u(t−1)−_ρs[1−_εss(5)−_η0(5)u(t)−_η1(5)u(t−1)]^{[_us(5)+_γ0(5)u(t)+_γ1(5)u(t−1)]2}−g_ρg∫05(_εss(y)+_η0(y)u(t)+_η1(y)u(t−1))dy−g_ρs∫05(1−_εss(y)−_η0(y)u(t)−_η1(y)u(t−1))dy.

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:

Δ_P0=_ρgUg2εss(0)+_η0(0)u(t)+_η1(0)u(t)+_ρs[1−_εss(0)−_η0(0)u(t)−_η1(0)u(t)]^{[_us(0)+_γ0(0)u(t)+_γ1(0)u(t)]2}−_ρgUg2εss(5)+_η0(5)u(t)+_η1(5)u(t)−_ρs[1−_εss(5)−_η0(5)u(t)−_η1(5)u(t)]^{[_us(5)+_γ0(5)u(t)+_γ1(5)u(t)]2}−g_ρg∫05(_εss(y)+_η0(y)u(t)+_η1(y)u(t))dy−g_ρs∫05(1−_εss(y)−_η0(y)u(t)−_η1(y)u(t))dy.

This is the pressure drop predicted by this model for constant inputu(t), 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ΔP(t)and the eventual steady pressure dropΔ_P0(t)is then equal to:

ΔP−Δ_P0=P(5)−P(0)−Δ_P0.

Linearizing Equation (72) aboutu(t)=u(t−1)gives:

ΔP−Δ_P0≈f(u(t−1))(u(t)−u(t−1)),

where:

f(y)=_k1(_k2+_k3y)^(_k2+_k4y)4+_k5 ^(_k6+_k7y)2+_k8(1−_k2−_k4y)(_k6+_k9y)+_k10(_k11+_k12y)^{(_k11+_k13y)4}+_k14 ^{(_k15+_k16y)2}+_k17(1−_k11−_k13y)(_k15+_k18y)+_k19,

and:

_k1=_{ρg Ug η1}(0), _k2=_εss(0), _k3=_η0(0)+2_η1(0),_k4=_η0(0)+_η1(0), _k5=_{ρs η1}(0), _k6=_us,ss(0),_k7=_γ0(0)+_γ1(0), _k8=−_{ρs γ1}(0), _k9=_γ0(0)+2_γ1(0),_k10=−_{ρg Ug η1}(5), _k11=_εss(5), _k12=_η0(5)+2_η1(5),_k13=_η0(5)+_η1(5), _k14=−_{ρs η1}(5), _k15=_us,ss(5),_k16=_γ0(5)+_γ1(5), _k17=_{ρs γ1}(5), _k18=_γ0(5)+2_γ1(5),_k19=g(_ρg−_ρs)_∫05η1(y)dy.

The input to the computational model was in terms of the velocity boundary condition, so thatu(t)=Δ_us(0,t). This can be connected to the inputs_S1and_S2via a quadratic model fitted to the test data where:

u(t)≈1_ε0(2_{a1 S1}(t−1)+_{a3 S2}(t)+_a4)(_S1(t)−_S1(t−1))+1_ε0(_{a1 S12}(t−1)+_{a2 S22}(t)+_{a3 S1}(t−1)_S2(t)+_{a4 S1}(t−1)+_{a5 S2}(t)+_a6)−_us,ss(0)=1_ε0(_{a1 S1}(t−1)+_{a3 S2}(t)+_a4)_S1(t)+1_ε0(_{a2 S22}(t)+_{a5 S2}(t)+_a6)−_us,ss(0).

Then:

ΔP−Δ_P0≈1_ε0f(u(t−1))(2_{a1 S1}(t−1)+_{a3 S2}(t)+_a4)_S1(t)+f(u(t−1))[1_ε0(_{a2 S22}(t)+_{a5 S2}(t)+_a6)−_us,ss(0)−u(t−1)]=_gΔP(_S1(t−1),_S2(t−1),_S2(t))_S1(t)+_fΔP(_S1(t−1),_S2(t−1),_S2(t)).

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

_yw(t)=f(y(t−1),⋯,y(t−_ny),_{S1_w}(t−1),⋯,_{S1_w}(t−_nu)),

where_S1,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)≈_yw(t)−_fΔP−_{gΔP S1}(t).

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:

_S1,fast(t)=_yr(t)−_yw(t)+_fΔPgΔP,

where_yr(t)is the reference signal. Hence, the final spatiotemporal wavelet controller_S1(t)implemented on the real CL process is given by:

_S1(t)=_S1,w(t)+_S1,fast(t).

This controller was also implemented in the single loop cold flow CL test rig._S1was taken as the single input andDP47as the output, while_S2 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∞ 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.

Figure 1. Alstom's combustion-gasification process.

Figure 2. Block diagram for a single-loop cold flow CL test rig.

Figure 3. Schemetic of the wavelet MRA-based self-tuning GPC control system.

Figure 4. Experimental facility for control testing.

Figure 5. PRBS test-input S1 and S2 (scfh).

Figure 6. PRBS test-Output DP47 (inch H2O).

Figure 7. PRBS test-Simulation data vs. experimental data.

Figure 8. PRBS test-One-step-ahead predictions vs. experimental data.

Figure 9. Pressure difference response of riser (DP47) during step setpoint changes.

Figure 10. Fluidizing air flow control (S1 and S2) during step setpoint changes.

Figure 11. Pressure difference response of riser during sinusoidal setpoint changes.

Figure 12. Fluidizing air flow control (S1 and S2) during setpoint changes.

Figure 13. Block diagram of controller implementation with fast dynamics.

Figure 14. Simulated impulse response of the 2-PDE riser model: (a)-the solid velocity response, (b)-the voidage response.

Figure 15. Wavelet-approximated impulse responseh(x,t): (a)-the solid velocity response approximation, (b)-the voidage response approximation.

Figure 16. Pressure difference response of riser (DP47) during step setpoint changes.

Figure 17. Fluidizing air flow control (S1 and S2) during setpoint changes.

Author Contributions

Conceptualization, S.Z., J.B., X.L. and C.N.; methodology, S.Z. and J.B.; software, S.Z. and X.L.; validation, S.Z., J.B., X.L., and C.N.; formal analysis, S.Z., J.B., Y.L., and H.E.-K.; investigation, S.Z., J.B., X.L., and C.N.; resources, J.B., X.L., and C.N.; data curation, S.Z. and X.L.; writing-original draft preparation, S.Z.; writing-review and editing, S.Z., J.B., X.L., C.N., Y.L., H.E.-K.; visualization, S.Z., J.B., X.L., and H.E.-K.; supervision, J.B., X.L., and C.N.; project administration, J.B., X.L., and C.N.; funding acquisition, C.N. All authors have read and agreed to the published version of the manuscript.

Funding

This work was funded by DOE/ALSTOM, contract number DE-FC26-07NT43095 and, in part, by by the National Institute of Biomedical Imaging and Bioengineering of the National Institutes of Health under award number R01EB029766. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.

Acknowledgments

The authors gratefully acknowledge the support of ALSTOM Power Plant Labs and DOE NETL. In particular, an invaluable organizational support of Ray P. Chamberland of ALSTOM and programmatic support of Robert R. Romanosky of DOE NETL is greatly appreciated. The authors are also greatly indebted to the DOE NETL project manager Susan M. Maley for her significant help with assessment of the state-of-the-art in mathematical modeling and numerical simulation of controlled CFBRs (continuous fluidized bed reactors). Cyrus Taft is gratefully acknowledged for facilitating participation of the University of Illinois team in the project. The University of Illinois team members Vivek Natarajan, currently at IIT Bombay, Bryan Petrus, currently at Nucor Steel, and Dong Ye, currently at Microsoft, are gratefully acknowledged for making innumerable contributions to all aspects of the project.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

References

1. Karthikeyan, A.; Previsic, M.; Scruggs, J.; Chertok, A. Non-linear Model Predictive Control of Wave Energy Converters with Realistic Power Take-off Configurations and Loss Model. In Proceedings of the 3rd IEEE Conference on Control Technology and Applications (CCTA 2019), Hong Kong, China, 19-21 August 2019; pp. 270-277.

2. Lou, X.; Neuschaefer, C.; Lei, H. Simulation and Advanced Controls for Alstom's Chemical Looping Process. In Proceedings of the 33rd International Technical Conference on Coal Utilization and Fuel Systems, Clearwater, FL, USA, 1-5 June 2008.

3. Lou, X.; Neuschaefer, C.; Lei, H. Simulation and Advanced Controls for Hybrid Combustion-Gasification Chemical Looping Process. In Proceedings of the 51st ISA Power Industry Division Symposium and 18th Annual Joint ISA POWID/EPRI Controls and Instrumentation Conference, Scottsdale, Arizona, 8-13 June 2008.

4. Fan, L.-S. Chemical Looping Systems for Fossil Energy Conversions; Wiley: Hoboken, NJ, USA, 2010.

5. Iloeje, C.O.; Zhao, Z.; Ghoniem, A.F. Design and techno-economic optimization of a rotary chemical looping combustion power plant with CO2 capture. Appl. Energy 2018, 231, 1179-1190.

6. Chen, P.; Sun, X.; Gao, M.; Ma, J.; Guo, Q. Transformation and migration of cadmium during chemical-looping combustion/gasification of municipal solid waste. Chem. Eng. J. 2019, 365, 389-399.

7. Vidal, R.; Bruna, J.; Giryes, R.; Soatto, S. Mathematics of Deep Learning. arXiv 2017, arXiv:1712.04741v1.

8. Bruna, J.; Mallat, S. Invariant Scattering Convolution Networks. IEEE Trans. Pattern Anal. Mach. Intell. 2013, 35, 1872-1886.

9. Alexandridis, A.K.; Zapranis, A.D. Wavelet neural networks: A practical guide. Neural Netw. 2013, 42, 1-27.

10. Doucoure, B.; Agbossou, K.; Cardenas, A. Time series prediction using artificial wavelet neural network and multi-resolution analysis: Application to wind speed data. Renew. Energy 2016, 92, 202-211.

11. Zhang, Q.; Benveniste, A. Wavelet networks. IEEE Trans. Neural Netw. 1992, 3, 889-898.

12. Zhang, Q. Using wavelet network in nonparametric estimation. IEEE Trans. Neural Netw. 1997, 8, 227-236.

13. Billings, S.A.; Wei, H.L. A new class of wavelet networks for nonlinear system identification. IEEE Trans. Neural Netw. 2005, 16, 862-874.

14. Coca, D.; Billings, A. Non-linear system identification using wavelet multiresolution models. Int. J. Control 2001, 74, 1718-1736.

15. Khelifa, S.; Gourine, B.; Rami, A.; Taibi, H. Assessment of nonlinear trends and seasonal variations in global sea level using singular spectrum analysis and wavelet multiresolution analysis. Arab. J. Geosci. 2016, 9, 560.

16. Sudre, J.; Yahia, H.; Pont, O.; Garçon, V. Ocean Turbulent Dynamics at Superresolution from Optimal Multiresolution Analysis and Multiplicative. IEEE Trans. Geosci. Remote Sens. 2015, 53, 6274-6285.

17. Hsu, C.; Lin, C.; Lee, T. Wavelet adaptive backstepping control for a class of nonlinear systems. IEEE Trans. Neural Netw. 2006, 17, 1175-1183.

18. Xu, J.; Tan, Y. Nonlinear adaptive wavelet control using constructive wavelet networks. IEEE Trans. Neural Netw. 2007, 18, 115-127.

19. Summers, S.; Jones, C.N.; Lygeros, J.; Morari, M. A Multiresolution Approximation Method for Fast Explicit Model Predictive Control. IEEE Trans. Autom. Control 2011, 56, 2530-2541.

20. Li, S.; Xi, Y. Application of wavelet to constrained generalized predictive control. In Proceedings of the 15th IEEE international Symposium on Intelligent Control (ISIC 2000), Rio, Patras, Greece, 17-19 July 2000; pp. 247-252.

21. Bitmead, R.; Gevers, M.; Wertz, V. Adaptive Optimal Control: The Thinking Man's GPC; Prentice Hall of Australia: Brunswick, Victoria, Australia, 1990.

22. Zhang, S.; Bentsman, J.; Lou, X.; Neuschaefer, C. Wavelet Multiresolution Model based Generalized Predictive Control for Hybrid Combustion-Gasification Chemical Looping Process. In Proceedings of the 51st IEEE Conference on Decision and Control, Maui, Hawaii, USA, 10-13 December 2012; pp. 2409-2414.

23. Bentsman, J.; Zhang, S.; Natarajan, V.; Petrus, B.; Ye, D. Development of Computational Approaches for Simulation and Advanced Controls for Hybrid Combustion-Gasification Chemical Looping; Technical Report; Department of Energy National Energy Technology Laboratory: Pittsburgh, PA, USA, 2011.

24. Zhang, S. Wavelet Adaptive and Predictive Control with Applications to Chemical Looping System. Ph.D. Thesis, University of Illinois at Urbana-Champaign, Urbana, IL, USA, 2014.

25. Clark, D.; Mohtadi, C.; Tuffs, P. Generalized predictive control-Part I. the basic algorithm. Automatica 1987, 23, 137-148.

26. Clark, D.; Mohtadi, C.; Tuffs, P. Generalized predictive control-Part II. Extensions and interpretations. Automatica 1987, 23, 149-160.

27. Matko, D.; Biasizzo, K.K. Generalized predictive control of a thermal plant using fuzzy model. In Proceedings of the 2000 American Control Conference, Chicago, IL, USA, 28-30 June 2000; Volume 3, pp. 2053-2057.

28. Gu, D.; Hu, H. Wavelet neural network-based predictive control for mobile robots. In Proceedings of the 2000 IEEE International Conference on Systems, Man and Cybernetics. 'Cybernetics Evolving to Systems, Humans, Organizations, and Their Complex Interactions', Nashville, TN, USA, 8-11 October 2000; Volume 5, pp. 3544-3549.

29. Yoo, S.J.; Choi, Y.H.; Park, J.B. Generalized predictive control based on self-recurrent wavelet neural network for stable path tracking of mobile robots: Adaptive learning rates approach. IEEE Trans. Circuits Syst. 2006, 53, 1381-1394.

30. Tao, J.; Dehmer, M.; Xie, G.; Zhou, Q. A Generalized Predictive Control-Based Path Following Method for Parafoil Systems in Wind Environments. IEEE Access 2019, 7, 42586-42595.

31. Li, S.; Jiang, P.; Han, K. Generalized Predictive Control of an Intensified Continuous Reactor. In Proceedings of the 37th Chinese Control Coference, Wuhan, China, 25-27 July 2018.

32. Grimble, M.J. Generalized predictive optimal control: An introduction to the advantages and limitations. Intern. J. Syst. Sci. 1992, 23, 85-98.

33. Mallat, J. A Theory for Multiresolution Signal Decomposition: The Wavelet Representation. IEEE Trans. Pattern Anal. Mach. Intell. 1989, 11, 674-693.

34. Chen, S.; Billings, S.A.; Luo, W. Orthogonal least squares methods and their application to nonlinear system identification. Int. J. Control 1989, 50, 1873-1896.

35. Astrom, K.J.; Wittenmark, B. Adaptive Control, 2nd ed.; Addison-Wesley: Reading, MA, USA, 1995.

36. Tsang, T.; Clarke, D. Generalized predictive control with input constraints. IEE Proc. 1988, 135, 451-460.

37. Huang, Y.; Turtona, R.; Park, J.; Famourib, P.; Boylec, E.J. Dynamic model of the riser in circulating fluidized bed. Powder Technol. 2006, 163, 23-31.

38. Narendra, K.S.; Mukhopadhyay, S. Adaptive control using neural networks and approximate models. IEEE Trans. Neural Netw. 1997, 8, 475-485.

39. Zhao, H.; Bentsman, J. Biorthogonal wavelet based identification of fast linear time varying systems-part 1: System representations. J. Dyn. Syst. Meas. Control 2001, 123, 585-592.

40. Zhao, H.; Bentsman, J. Biorthogonal wavelet based identification of fast linear time varying systems-part 2: Algorithms and performance analysis. J. Dyn. Syst. Meas. Control 2001, 123, 593-600.

41. Opfermann, J.D.; Leonard, S.; Decker, R.S.; Uebele, N.A.; Bayne, C.E.; Joshi, A.S.; Krieger, A. Semi-Autonomous Electrosurgery for Tumor Resection Using a Multi-Degree of Freedom Electrosurgical Tool and Visual Servoin. In Proceedings of the 2017 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), Vancouver, BC, Canada, 24-28 September 2017; pp. 3653-3660.

42. Madhukar, A.; Park, Y.; Kim, W.; Sunaryanto, H.J.; Berlin, R.; Chamorro, L.P.; Bentsman, J.; Ostoja-Starzewski, M. Heat conduction in porcine muscle and blood: Experiments and time-fractional telegraph equation model. J. R. Soc. Interface 2019, 16, 1-8.