An Active Fault-Tolerant PWM Tracker for Unknown Nonlinear Stochastic Hybrid Systems: NARMAX Model and OKID-Based State-Space Self-Tuning Control

Abstract

An active fault-tolerant pulse-width-modulated tracker using the nonlinear autoregressive moving average with exogenous inputs model-based state-space self-tuning control is proposed for continuous-time multivariable nonlinear stochastic systems with unknown system parameters, plant noises, measurement noises, and inaccessible system states. Through observer/Kalman filter identification method, a good initial guess of the unknown parameters of the chosen model is obtained so as to reduce the identification process time and enhance the system performances. Besides, by modifying the conventional self-tuning control, a fault-tolerant control scheme is also developed. For the detection of fault occurrence, a quantitative criterion is exploited by comparing the innovation process errors estimated by the Kalman filter estimation algorithm. In addition, the weighting matrix resetting technique is presented by adjusting and resetting the covariance matrix of parameter estimates to improve the parameter estimation for faulty system recovery. The technique can effectively cope with partially abrupt and/or gradual system faults and/or input failures with fault detection.

1. Introduction

One of the major challenges in the identification of nonlinear stochastic systems is the choice of model representation and structure. Generally speaking, polynomial expressions are extensively used to represent nonlinear dynamic systems when the response of a system is dominated by nonlinear characteristics, so it is usually necessary to use an adequate nonlinear model to analyze and synthesize nonlinear systems. The Nonlinear AutoRegressive Moving Average with eXogenous inputs (NARMAX) model was first introduced and rigorously derived by [1, 2], which provides a unified representation for a wide class of nonlinear stochastic systems. The initial parameter of NARMAX model is important to reduce the time of the identifying process; therefore observer/Kalman filter identification (OKID) [3] is applied to estimate the initial parameter of NARMAX model in this paper.

Over the past decades, the state-space self-tuning control (STC) methods initially introduced in [4, 5] have been shown to be effective in designing advanced adaptive controllers for linear multivariable stochastic systems [6]. In those approaches, the illustrious Kalman state-estimation algorithm [7] has been embedded into an online parameter estimation algorithm. They utilize state-space self-tuners based on innovation models, where (i) the equivalent internal states can be estimated successively; (ii) the stable/unstable and minimum/nonminimum-phase multivariable systems can be controlled accurately; (iii) the self-tuners are simple, reliable, and robust; (iv) the adaptive Kalman gain can subsequently be obtained.

In the state-space self-tuning control of a class of multivariable stochastic systems, it is often required to transform the state equation in an observer block companion form to the equivalent state equation in a controller block companion form. However, this method involves the products of high-dimensional system matrices; hence, it may introduce large computational errors and should be avoided. In this paper, the suboptimal tracker is directly designed in an observer block companion form. One point must be noticed that the state-space self-tuning control scheme for nonlinear stochastic hybrid systems proposed by Guo et al. [8] estimates the system parameters at every sampling instant, and then an adaptive controller could be designed based on the parameter estimated at every sampling instant. The framework of the state-space STC seems to agree with that of the active fault tolerance in a real time. For faulty system recovery, the modified Kalman filter estimation algorithm can be used to estimate system parameters, instead of using the recursive extended-least-squares (RELSs) algorithm in the conventional STC scheme. An important procedure in the proposed fault-tolerant control (FTC) is to determine the initializations of the covariance matrices of process noises, measurement noises, and parameter estimation errors and to reset these matrices for improvement of the parameter estimation [9] whenever the unanticipated faults happen. The fault detection is decided by the given preset threshold. As for the faults, both abrupt and gradual faults are considered in this paper.

This paper is organized as follows. In Section 2, the self-tuning control of stochastic systems based on NARMAX model is introduced. Section 3 proposes an online identification and observer with OKID for self-tuning control. Section 4 discusses the NARMAX model-based state-space self-tuner for unknown nonlinear stochastic hybrid systems. In Section 5, a fault-tolerant scheme is proposed by modifying the conventional state-space self-tuning control approach for unknown multivariable stochastic systems. Finally, two illustrative examples are shown in Section 6.

2. Self-Tuning Control Scheme Based on NARMAX Model

2.1. Digital Linear Suboptimal Tracker Design

Consider the underlying linear discrete-time system:𝑥((𝑘+1)𝑇)=𝐺𝑥(𝑘𝑇)+𝐻𝑢(𝑘𝑇),𝑦(𝑘𝑇)=𝐶𝑥(𝑘𝑇), (1a)(1b)where 𝑥(𝑘𝑇)∈𝑅^𝑛, 𝑢(𝑘𝑇)∈𝑅^𝑚, and 𝑦(𝑘𝑇)∈𝑅^𝑝 are the state variable, control input, and control output, respectively. 𝑇 stands for the sampling period. For simplifying equations and causing no ambiguity, •(𝑘𝑇) will be replaced by •(𝑘) in many equations of this paper.

It is desired to minimize the following quadratic cost function:𝐽_𝑑=12{^∞_𝑘=1[𝐶𝑥(𝑘)−Γ(𝑘)]^𝑇𝑄_𝑑[𝐶𝑥(𝑘)−Γ(𝑘)]+𝑢^𝑇(𝑘)𝑅_𝑑𝑢(𝑘)}, (2) where 𝑄_𝑑≥0, 𝑅_𝑑>0 and Γ(𝑘)∈𝑅^𝑝 is a reference input. The superscript 𝑇 designates the transpose of a matrix or vector. For the linear discrete-time system (1a) and (1b), the digital linear suboptimal tracker design is given by𝑢(𝑘)=−𝐾_𝑑𝑥(𝑘)+𝐸_𝑑Γ^∗(𝑘), (3) where𝐾_𝑑=(𝐻^𝑇𝑃𝐻+𝑅_𝑑)⁻¹𝐻^𝑇𝐸𝑃𝐺,_𝑑=(𝐻^𝑇𝑃𝐻+𝑅_𝑑)⁻¹𝐻^𝑇[𝐼−(𝐺−𝐻𝐾_𝑑)^𝑇]⁻¹𝐶^𝑇𝑄_𝑑. (4)Γ^∗(𝑘) = Γ(𝑘+1) for the tracking purpose, and 𝑃 is the positive definite and symmetric solution to the following Riccati equation:𝐺^𝑇𝑃𝐺−𝑃−(𝐺^𝑇𝑃𝐻)(𝐻^𝑇𝑃𝐻+𝑅_𝑑)⁻¹(𝐻^𝑇𝑃𝐺)+𝐶^𝑇𝑄_𝑑𝐶=0. (5) If the discrete-time system is time-varying, the digital tracker gains (𝐾_𝑑,𝐸_𝑑) would also be time-varying.

2.2. Design of PWM Controller from PAM Controller

In general, there exist two types of digital controllers: the pulse-amplitude-modulated (PAM) controller and the pulse-width-modulated (PWM) controller [10, 11]. The PAM controller, which produces a series of piecewise-constant continuous pulses having variable amplitude and fixed width, is commonly utilized in digital control of all types. The PWM controller, which produces a series of discontinuous pulses with a fixed amplitude and variable width, has become popular in industry for on-off control of DC power converters and stepper motors (widely used in robotics), satellite station-keeping (with on-off reaction jets), and so forth. Since the conventional direct digital design approach takes into account only the sampling instants of the continuous-time system, the resulting PAM and PWM controllers could produce degradation in the intersample behavior of the closed-loop sampled-data system. Hence, we focus on the digital approach for the development of the PAM and PWM controllers.

The PAM controlled sampled-data system is described by𝑥̇(𝑡)=𝐴𝑥(𝑡)+^𝑚_𝑗=1𝐵^(𝑗)𝑢^(𝑗)(𝑘𝑇), (6) where 𝐵^(𝑗) is the 𝑗th column of 𝐵 and 𝑢^(𝑗)(𝑘𝑇) is the 𝑗th component of the PAM input vector. The corresponding discrete-time model of the system (6) with a zero-order hold is𝑥(𝑘+1)=𝐺𝑥(𝑘)+^𝑚_𝑗=1𝐻^(𝑗)𝑢^(𝑗)(𝑘), (7a) where𝐺=𝑒^𝐴𝑇,𝐻^(𝑗)=(𝐺−𝐼)𝐴⁻¹𝐵^(𝑗),if𝐴⁻¹exists,_∞𝑖=0𝑇^𝑖+1𝐴(𝑖+1)!^𝑖𝐵^(𝑗),if𝐴⁻¹doesnotexist. (7b)Here 𝐼 denotes the identity matrix with the same size as 𝐺.

Let a PWM controlled sampled-data system be represented as𝑥̇_𝑤(𝑡)=𝐴𝑥_𝑤(𝑡)+^𝑚_𝑗=1𝐵^(𝑗)𝑢_𝑤^(𝑗)(𝜆), (8) where the PWM controller is𝑢_𝑤^(𝑗)(𝜆)=0,𝑘𝑇≤𝜆<𝑘𝑇+𝜏_𝑘^(𝑗),𝑢_𝑀^(𝑗),𝑘𝑇+𝜏_𝑘^(𝑗)≤𝜆<𝑘𝑇+𝜏_𝑘^(𝑗)+𝛿_𝑘^(𝑗),0,𝑘𝑇+𝜏_𝑘^(𝑗)+𝛿_𝑘^(𝑗)≤𝜆<𝑘𝑇+𝑇,for𝑗=1,2,…,𝑚. (9) Here 𝑢_𝑀^(𝑗), 𝜏_𝑘^(𝑗), and 𝛿_𝑘^(𝑗) are the amplitude of the 𝑗th component of the input vector, the firing delay, and the firing duration in the PWM mode at the 𝑘th time step, respectively. The graphical illustration of PAM and PWM inputs is shown in Figure 1. The values of 𝜏_𝑘^(𝑗) and 𝛿_𝑘^(𝑗) can be determined as follows:𝛿_𝑘^(𝑗)=𝑇𝑢^(𝑗)(𝑘𝑇)𝑢_𝑀^(𝑗),𝜏_𝑘^(𝑗)=(𝑇−𝛿_𝑘^(𝑗))2. (10) Therefore, the aforementioned PAM controller (3) can be equivalently transferred into the corresponding PWM controller (9) by (10).

Figure 1 

Graphical illustration of PAM and equivalent PWM inputs.

2.3. The Structure of the State-Space STC

The structure of the state-space STC scheme includes a parameter and state estimator and a controller design. A typical state-space STC structure is illustrated in Figure 2.

Figure 2 

Block diagram of a typical state-space self-tuning control.

Under this framework, the parameters and states of the unknown system are estimated from the consecutive control input data (𝑢(𝑘),𝑢(𝑘−1),𝑢(𝑘−2),…) and the system output data (𝑦(𝑘),𝑦(𝑘−1),𝑦(𝑘−2),…). To be useful in STC the parameter estimation scheme should be iterative, allowing the estimated model of the controlled system to be updated at each sample interval until the control goal is achieved. In this scheme new input/output data become available at each sample interval. First, assume that the parameter estimate 𝜃(𝑘−1) is obtained based on the past information up to the time step 𝑘−1. This could be used to yield an estimate of the output at the current time step 𝑘. The estimate of the current output is then compared with the observed current output to generate a prediction error. This in turn generates an update to the model which corrects 𝜃(𝑘−1) to the new value 𝜃(𝑘). Depending on the parameter estimate 𝜃(𝑘), appropriate controllers can be designed. The reference input Γ(𝑘) and the designed adaptive controller generate real-time control actions for unknown dynamic systems. Notice that if the control input is persistently excited, then the convergence to the true system parameters is guaranteed [12].

2.4. NARMAX Model for Self-Tuning Control Scheme

The ARMAX model-based state-space self-tuning control has been represented in [8]. Nevertheless, it is well known that the NARMAX model is a general and natural representation of nonlinear systems; however, the NARMAX model-based state-space self-tuner for fault-tolerant control has not been developed in literature. Naturally, the objective of this paper is to extend the ARMAX model-based state-space self-tuner to the NARMAX model-based state-space self-tuner for a high-performance tracking control.

The expression of the NARMAX model proposed here for the m-input p-output system is given by𝑦_𝑖(𝑘)=𝐹_𝑖(𝑦(𝑘−1),…,𝑦(𝑘−𝑛_𝑦),𝑢(𝑘−1),…,𝑢(𝑘−𝑛_𝑢),𝜀(𝑘−1),…,𝜀(𝑘−𝑛_𝜀))+𝜀_𝑖=(𝑘)_{𝜉𝑖𝑗=1}𝜃_𝑖𝑗𝜙_𝑖𝑗(𝑘)+𝜀_𝑖(𝑘)=𝜃^𝑇_𝑖𝜙_𝑖(𝑘)+𝜀_𝑖(𝑘),for𝑖=1,2,…,𝑝. (11) Here 𝑦_𝑖(𝑘)∈ℜ is the 𝑖th output at time step 𝑘 (𝑘=1,2,3,…); 𝑢(𝑘)∈ℜ^𝑚 is the control input; notation 𝜀_𝑖(𝑘)∈ℜ denotes the ith component of the prediction error 𝜀(𝑘)∈ℜ^𝑝. When the true system output is corrupted by white or colored noise, knowledge of the past values of the modeling error is required to estimate unknown parameters. However, the modeling error is, in general, an unobservable error process, since it depends on the unknown noise. A number of estimation procedures exist, nevertheless, which replace the modeling error by an estimate, usually taken to be the prediction error or the residual. In this paper the RELS algorithm will be utilized to estimate unknown parameters of the NARMAX model; hence, the prediction error 𝜀(𝑘) will be instead of the modeling error. In addition, 𝑛_𝑦, 𝑛_𝑢, and 𝑛_𝜀 are the maximum lag indices of 𝑦(𝑘), 𝑢(𝑘), and 𝜀(𝑘), respectively. 𝜃_𝑖 is the unknown parameter vector for the 𝑖th output. 𝜉_𝑖 denotes the number of 𝜃_𝑖. The vector 𝜙_𝑖(𝑘) is the linear or nonlinear function of (𝑦(𝑘−1),…,𝑦(𝑘−𝑛_𝑦),𝑢(𝑘−1),…,𝑢(𝑘−𝑛_𝑢),𝜀(𝑘−1),…,𝜀(𝑘−𝑛_𝜀)) related to the ith output. Although 𝐹_𝑖(⋅) may be linear or nonlinear polynomial in terms of the past values of (𝑦(𝑘),𝑢(𝑘),𝜀(𝑘)), the whole model is nonlinear. 𝜙_𝑖(𝑘) may include any linear or nonlinear variables such as terms raised to an integer of power (e.g., 𝑦³₁(𝑘−5), 𝑢²₃(𝑘−4)), products of present and past inputs (e.g., 𝑢(𝑘)𝑢(𝑘−2)), products of past outputs (e.g., 𝑦₂(𝑘−2)𝑦₃(𝑘−4)), or cross-terms (e.g., 𝑦₄(𝑘−1)𝑢₂(𝑘−3)).

To determine the unknown parameter vector 𝜃_𝑖, the standard RELS algorithm can be applied to the linear-in-the-parameter model (11). This RELS algorithm is given by𝜆(𝑘)=𝜆₀𝜆(𝑘−1)+(1−𝜆₀𝜀),_𝑖(𝑘)=𝑦_𝑖(𝑘)−𝜃^𝑇_𝑖(𝑘−1)𝜙_𝑖𝜃(𝑘),_𝑖(𝑘)=𝜃_𝑖𝑆(𝑘−1)+^𝑖(𝑘−1)𝜙_𝑖(𝑘)𝜆(𝑘)+𝜙_𝑖^𝑇(𝑘)𝑆_𝑖(𝑘−1)𝜙_𝑖𝜀(𝑘)^𝑖𝑆(𝑘),_𝑖1(𝑘)=𝜆(𝑘){𝑆_𝑖𝑆(𝑘−1)−^𝑖(𝑘−1)𝜙_𝑖(𝑘)𝜙^𝑇_𝑖(𝑘)𝑆_𝑖(𝑘−1)𝜆(𝑘)+𝜙^𝑇_𝑖(𝑘)𝑆_𝑖(𝑘−1)𝜙_𝑖(𝑘)}. (12)(1)(2)(3) Here 𝜆(𝑘) is referred to as the forgetting factor to discount the old measurements with the initial condition 0.9<𝜆(0)<1 and the updating factor 0<𝜆₀<1, while 𝜆(𝑘)=1 recovers the original RELS algorithm without forgetting where all data were weighted equally no matter how far back in the past. The adjustment of 𝜆(𝑘) is a tradeoff between high robustness against noises (large 𝜆(𝑘)) and fast tracking capability (small 𝜆(𝑘)). The prediction error 𝜀_𝑖(𝑘) is the difference between the measured output 𝑦_𝑖(𝑘) and the one-step-ahead prediction 𝑦_𝑖(𝑘∣𝑘−1)=𝜃^𝑇_𝑖(𝑘−1)𝜙_𝑖(𝑘) of 𝑦_𝑖(𝑘) made at time step 𝑘−1 based on the model corresponding to the estimate 𝜃_𝑖(𝑘−1).

If the prediction error is small, the estimate 𝜃_𝑖(𝑘−1) is “good” and should not be modified very much. The term 𝑆_𝑖(𝑘) is the parameter estimation error covariance matrix with 𝑆_𝑖(0)=𝛾_𝑖𝐼_𝜉𝑖, where 𝐼_𝜉𝑖 is called the identity matrix of order 𝜉_𝑖, and 𝛾_𝑖 is a sufficiently large positive number in order that the parameter estimate can quickly jump away from 𝜃_𝑖(0).

Owing to various kinds of combinations of (𝑦(𝑘),𝑢(𝑘),𝜀(𝑘)) for 𝐹_𝑖(⋅), many classes of NARMAX models can be chosen. To simplify the whole control scheme, it is desired to choose some simple structures of dynamic nonlinear models. Thus, the state-space STC scheme with the NARMAX model for nonlinear stochastic systems can work fast and more precisely.

3. The Fast Online Identification and Observer Design Based on OKID

In this section, we slightly modify the basic structure of discrete-time state-space model, which is useful for the hybrid state-space self-tuning control law design for both linear and nonlinear systems. Here we consider a class of multivariable nonlinear systems. Once having the parameter estimate 𝜃_𝑖(𝑘) from the standard RELS algorithm, the NARMAX model for the STC scheme can accurately approximate the responses of the nonlinear system. Moreover, the initial parameter 𝜃_𝑖(0) of the NARMAX model will affect the convergent speed of RELS process. In order to get a suitable initial parameter 𝜃_𝑖(0) to shorten the transient process of RELS, one could apply the observer/Kalman filter identification (OKID) to evaluate it here.

The regression vector 𝜙_𝑖(𝑘) in (11) is composed of𝑦^𝑇𝒜𝑦^𝑇ℬ⋯𝑦^𝑇(𝑘−𝑛_𝑦)𝑢^𝑇𝒜𝑢^𝑇ℬ⋯𝑢^𝑇(𝑘−𝑛_𝑢)𝜀^𝑇𝒜𝜀^𝑇ℬ⋯𝜀^𝑇(𝑘−𝑛_𝜀), (13)where 𝒜 denotes (𝑘−1) and ℬ denotes (𝑘−2). They are not independent factors, so it is difficult to design a digital controller directly from the STC scheme with the NARMAX model. For this reason, one could apply the optimal linearization to the NARMAX model to configure a linear discrete-time state-space observer so as to design a digital linear controller of the STC scheme. Thus the well-designed control law 𝑢(𝑘) could make the output 𝑦(𝑘) of the unknown system track exactly the prespecified reference Γ(𝑘).

3.1. Preliminary of Discrete State-Space Observer

A preliminary structure of the discrete state-space observer of the linear system is presented in [6]. Consider the linear discrete-time stochastic system characterized by𝑥(𝑘+1)=𝐺𝑥(𝑘)+𝐻𝑢(𝑘)+𝑤(𝑘),𝑦(𝑘)=𝐶𝑥(𝑘)+𝑣(𝑘), (14) where 𝐺∈ℜ^𝑛×𝑛, 𝐻∈ℜ^𝑛×𝑚, and 𝐶∈ℜ^𝑝×𝑛 are system, input, and output matrices, respectively; 𝑥(𝑘)∈ℜ^𝑛, 𝑢(𝑘)∈ℜ^𝑚, and 𝑦(𝑘)∈ℜ^𝑝 are state, input, and output vectors, respectively; 𝑤∈ℜ^𝑛 and 𝑣∈ℜ^𝑝 are assumed to be stationary white noise processes with zero mean values and the covariance matrix:𝐸(𝑤𝑤(𝑘)𝑣(𝑘)^𝑇(𝑙)𝑣^𝑇)=𝑆(𝑙)𝑄𝑆_𝑇𝑅𝛿(𝑘−𝑙),for𝑘,𝑙=0,1,2,…, (15) where 𝐸(⋅) is termed the expected value of a random sequence {⋅}. 𝑄 and 𝑅 are nonnegative and positive definite symmetric matrices for all 𝑘, respectively. Due to the assumption that both white noise processes are stationary, the covariance matrix of the process noise, 𝑄, and the covariance matrix of the measurement noise, 𝑅, are constant; otherwise, they might change with each time step or measurement. The symbol 𝛿(𝜏) is the Kronecker delta function, which is 1 for 𝜏=0 and 0 for 𝜏≠0. The system (14) is said to be in the block observable form, if the following observability matrixΘ=(𝐶𝐺^𝑟−1)^𝑇(𝐶𝐺^𝑟−2)^𝑇⋯(𝐶𝐺)^𝑇𝐶^𝑇𝑇 (16) is of full rank.

Note that the observability index of Θ is 𝑟=𝑛/𝑝, if it is an integer (otherwise, it is undefined). This constraint means that the Kronecker indices of the system are all such integer 𝑟 that satisfies 𝑛=𝑟𝑝. When this system (14) is block observable, it can be transformed into the block observable companion form as follows:𝑥_𝑜(𝑘+1)=𝐺_𝑜𝑥_𝑜(𝑘)+𝐻_𝑜𝑢(𝑘)+𝑤_𝑜(𝑘),𝑦(𝑘)=𝐶_𝑜𝑥_𝑜(𝑘)+𝑣_𝑜(𝑘), (17) where𝑇_𝑜1=Θ⁻¹𝐼_𝑝0_𝑝⋯0_𝑝^𝑇∈ℜ^𝑛×𝑝,𝑇_𝑜=𝐺^𝑟−1𝑇_𝑜1𝐺^𝑟−2𝑇_𝑜1⋯𝐺𝑇_𝑜1𝑇_𝑜1∈ℜ^𝑛×𝑛,𝑥_𝑜(𝑘)=𝑥^𝑇_𝑜1(𝑘)𝑥^𝑇_𝑜2(𝑘)⋯𝑥^𝑇_𝑜𝑟(𝑘)^𝑇=𝑇_𝑜⁻¹𝐺𝑥(𝑘),_𝑜=𝑇_𝑜⁻¹𝐺𝑇_𝑜=−𝐺_𝑜1𝐼_𝑝0_𝑝⋯0_𝑝−𝐺_𝑜20_𝑝𝐼_𝑝⋯0_𝑝⋮⋮⋮⋱⋮−𝐺_{𝑜(𝑟−1)}0_𝑝0_𝑝⋯𝐼_𝑝−𝐺_𝑜𝑟0_𝑝0_𝑝⋯0_𝑝,𝐻_𝑜=𝑇_𝑜⁻¹𝐻=𝐻^𝑇_𝑜1𝐻^𝑇_𝑜2⋯𝐻^𝑇_𝑜𝑟^𝑇,𝐶_𝑜=𝐶𝑇_𝑜=𝐶_𝑜1𝐶_𝑜2⋯𝐶_𝑜𝑟=𝐼_𝑝0_𝑝⋯0_𝑝∈ℜ^{𝑝×𝑟𝑝},𝑤_𝑜(𝑘)=𝑇_𝑜⁻¹𝑣𝑤(𝑘),_𝑜(𝑘)=𝑣(𝑘),𝑤𝐸{^𝑜𝑣(𝑘)_𝑜𝑤(𝑘)^𝑇𝑜(𝑙)𝑣^𝑇_𝑜}=𝑇(𝑙)_𝑜−1𝑄(𝑇_𝑜⁻¹)^𝑇𝑇_𝑜⁻¹𝑆(𝑇_𝑜⁻¹𝑆)^𝑇𝑅≜𝑄𝛿(𝑘−𝑙)_𝑜𝑆_𝑜𝑆^𝑇_𝑜𝑅_𝑜𝛿(𝑘−𝑙), (18) in which 0_𝑝 and 𝐼_𝑝 stand for the 𝑝×𝑝 null matrix and identity matrix, respectively.

The system (17) can be represented by a state-space innovation model [13] as follows:𝑥_𝑜(𝑘+1∣𝑘)=𝐺_𝑜𝑥_𝑜(𝑘∣𝑘−1)+𝐻_𝑜𝑢(𝑘)+𝐾_𝑜(𝑘)𝑒_𝑜𝑒(𝑘),_𝑜(𝑘)=𝑦(𝑘)−𝐶_𝑜𝑥_𝑜(𝑘∣𝑘−1). (19) Here the Kalman gain 𝐾_𝑜(𝑘) can be computed by the following algorithm [14]:𝐾_𝑜(𝑘)=[𝐺_𝑜𝐽_𝑜(𝑘)𝐶^𝑇_𝑜+𝑆_𝑜][𝐶_𝑜𝐽_𝑜(𝑘)𝐶^𝑇_𝑜+𝑅_𝑜]⁻¹,𝐽_𝑜(𝑘+1)=[𝐺_𝑜−𝐾_𝑜(𝑘)𝐶_𝑜]𝐽_𝑜(𝑘)[𝐺_𝑜−𝐾_𝑜(𝑘)𝐶_𝑜]^𝑇+𝐾_𝑜(𝑘)𝑅_𝑜𝐾^𝑇_𝑜(𝑘)−𝑆_𝑜𝐾^𝑇_𝑜(𝑘)−𝐾_𝑜(𝑘)𝑆^𝑇_𝑜+𝑄_𝑜=𝐸{[𝑥_𝑜(𝑘+1)−𝑥_𝑜(𝑘+1∣𝑘)]×[𝑥_𝑜(𝑘+1)−𝑥_𝑜(𝑘+1∣𝑘)]^𝑇𝐽},_𝑜(0)=𝐸{[𝑥_𝑜(0)−𝑥_𝑜(0)][𝑥_𝑜(0)−𝑥_𝑜(0)]^𝑇}, (20) in which 𝑥_𝑜(𝑘∣𝑘−1) is the optimal estimate of 𝑥_𝑜(𝑘) given by the measurement data 𝑦(𝑖) and control input data 𝑢(𝑖) for 𝑖=0,1,…,𝑘−1; 𝑒_𝑜(𝑘) is called the innovation process with zero mean and the covariance matrix 𝐸{𝑒_𝑜(𝑘)𝑒^𝑇_𝑜(𝑘)}=𝐶_𝑜𝐽_𝑜(𝑘)𝐶^𝑇_𝑜+𝑅_𝑜.

If the pair (𝐺_𝑜,𝐶_𝑜) is detectable and the pair (𝐺_𝑜−𝑆_𝑜𝑅_𝑜⁻¹𝐶_𝑜,𝑄_𝑜) is completely stabilizable for any 𝑄_𝑜 satisfying with 𝑄_𝑜𝑄^𝑇_𝑜=𝑄_𝑜−𝑆_𝑜𝑅_𝑜⁻¹𝑆^𝑇_𝑜 [15], then 𝐽_𝑜(𝑘)→𝐽_𝑜(∞), where 𝐽_𝑜(∞) is the stationary error covariance matrix, so that 𝐾_𝑜(𝑘)→𝐾_𝑜(∞) (the stationary Kalman gain) as 𝑘→∞. Furthermore, the eigenvalues of 𝐺_𝑜−𝐾_𝑜(∞)𝐶_𝑜 are all inside the unit circle centered at the origin. Let 𝑧⁻¹ be the backward time-shift operator and set 𝐾_𝑜=[𝐾^𝑇_𝑜1,𝐾^𝑇_𝑜2,…,𝐾^𝑇_𝑜𝑟]^𝑇∈ℜ^𝑛×𝑝. Then, the input-output relationship of the steady-state innovation representation of (19) can be rewritten as 𝑦(𝑘)=𝐶_𝑜[𝐼_𝑛−𝐺_𝑜𝑧⁻¹]⁻¹[𝐻_𝑜𝑧⁻¹𝑢(𝑘)+𝐾_𝑜(∞)𝑧⁻¹𝑒_𝑜(𝑘)]+𝑒_𝑜(𝑘)=[𝐺_𝑙(𝑧⁻¹)]⁻¹[𝐻_𝑙(𝑧⁻¹)]𝑢(𝑘)+[𝐺_𝑙(𝑧⁻¹)]⁻¹[𝐷_𝑙(𝑧⁻¹)]𝑒_𝑜(𝑘), (21a) or𝐺_𝑙(𝑧⁻¹)𝑦(𝑘)=𝐻_𝑙(𝑧⁻¹)𝑢(𝑘)+𝐷_𝑙(𝑧⁻¹)𝑒_𝑜(𝑘), (21b)where𝐺_𝑙(𝑧⁻¹)=𝐼_𝑝+𝐺_𝑜1𝑧⁻¹+𝐺_𝑜2𝑧⁻²+⋯+𝐺_𝑜𝑟𝑧^−𝑟,𝐻_𝑙(𝑧⁻¹)=𝐻_𝑜1𝑧⁻¹+𝐻_𝑜2𝑧⁻²+⋯+𝐻_𝑜𝑟𝑧^−𝑟,𝐷_𝑙(𝑧⁻¹)=𝐼_𝑝+𝐷_𝑜1𝑧⁻¹+𝐷_𝑜2𝑧⁻²+⋯+𝐷_𝑜𝑟𝑧^−𝑟,𝐷_𝑜𝑖=𝐺_𝑜𝑖+𝐾_𝑜𝑖,𝑖=1,2,…,𝑟. (22)

Note that all the zeros of det(𝐷_𝑙(𝑧⁻¹)) must be inside the unit circle centered at the origin. Notation det(⋅) means the determinant of ⋅. From (21b) we can easily see that (21a) is a multivariable ARMAX model. If these parameter matrices 𝐺_𝑜, 𝐻_𝑜, and 𝐶_𝑜 in (17) are known, and the covariance matrix therein is available, then the recursive estimation algorithm (20) can be applied to determine the Kalman gain 𝐾_𝑜(𝑘). Thus, the state 𝑥_𝑜(𝑘) can be optimally estimated by using the state-space innovation model (19). When the matrices 𝐺_𝑜 and 𝐻_𝑜 are unknown, and the covariance matrix is unavailable, the following model in the block observable form𝑥_𝑜(𝑘+1∣𝑘)=𝐺_𝑜(𝑘∣𝑘−1)𝑥_𝑜+𝐻(𝑘∣𝑘−1)_𝑜(𝑘)𝑢(𝑘)+𝐾_𝑜(𝑘)𝑒_𝑜(𝑘),𝑒_𝑜(𝑘)=𝑦(𝑘)−𝐶_𝑜𝑥_𝑜(𝑘∣𝑘−1), (23) can be used in conjunction with the RELS algorithm to find the Kalman gain estimate 𝐾_𝑜(𝑘) and the state estimate 𝑥_𝑜(𝑘). Here 𝐺_𝑜(𝑘), 𝐻_𝑜(𝑘), and 𝐾_𝑜(𝑘) contain parameter estimates 𝐺_𝑜𝑖(𝑘), 𝐻_𝑜𝑖(𝑘), and 𝐾_𝑜𝑖(𝑘) for 𝑖 = 1,2,…,𝑟. When all these parameter estimates converge to the true values, the state estimate 𝑥_𝑜(𝑘∣𝑘−1) of the state 𝑥_𝑜(𝑘) and the innovation process 𝑒_𝑜(𝑘) converge to the optimal state estimate 𝑥_𝑜(𝑘∣𝑘−1) and innovation process 𝑒_𝑜(𝑘), respectively. The linear regression model of (23) can be given in the following form:𝑦(𝑘)=𝜃^𝑇(𝑘−1)𝜙(𝑘)+𝜀(𝑘), (24a) where𝜃(𝑘)=𝐺_𝑜1(𝑘)𝐺_𝑜2(𝑘)⋯𝐺_𝑜𝑟(𝑘)𝐻_𝑜1(𝑘)𝐻_𝑜2(𝑘)⋯𝐻_𝑜𝑟(𝑘)𝐷_𝑜1(𝑘)𝐷_𝑜2(𝑘)⋯𝐷_𝑜𝑟(𝑘)^𝑇,𝜙(𝑘)=−𝑦^𝑇𝒜−𝑦^𝑇ℬ⋯−𝑦^𝑇(𝑘−𝑟)𝑢^𝑇𝒜𝑢^𝑇ℬ⋯𝑢^𝑇(𝑘−𝑟)𝜀^𝑇𝒜𝜀^𝑇ℬ⋯𝜀^𝑇(𝑘−𝑟)^𝑇. (24b)

It is important to note that as long as the matrix𝐺_𝑐=𝐺_𝑜−𝐾_𝑜𝐶_𝑜=−𝐷_𝑜1𝐼_𝑝0_𝑝⋯0_𝑝−𝐷_𝑜20_𝑝𝐼_𝑝⋯0_𝑝⋮⋮⋮⋱⋮−𝐷_{𝑜(𝑟−1)}0_𝑝0_𝑝⋯𝐼_𝑝−𝐷_𝑜𝑟0_𝑝0_𝑝⋯0_𝑝 (25) is asymptotically stable, the boundedness of the noise sequences implies that the estimation error will always be bounded. Whenever 𝐷_𝑜𝑖(𝑘)=0_𝑝, for 𝑖=1,2,…,𝑟, it designates a dead-beat-like property.

3.2. OKID Formulation

Consider again the linear discrete-time system represented by (1a) and (1b). When states of an observable system are inaccessible, an observer is usually applied to estimate the states from the information of input and output. Therefore, given the observer gain 𝐹, the system (1a) can be rewritten as=𝑥(𝑘+1)=𝐺𝑥(𝑘)+𝐻𝑢(𝑘)+𝐹𝑦(𝑘)−𝐹𝑦(𝑘)=(𝐺+𝐹𝐶)𝑥(𝑘)+𝐻𝑢(𝑘)−𝐹𝑦(𝑘)𝐺𝑥(𝑘)+𝐻𝑣(𝑘), (26a) where𝐺=𝐺+𝐹𝐶,𝐻=𝐻−𝐹,𝑣(𝑘)=𝑢(𝑘)𝑦(𝑘). (26b)(26c)(26d)

For a nonzero initial condition, 𝑥(0)≠0, the approach should be taken. For this purpose, (26a)–(26d) are extended into the following form:𝑥(𝑘+1)=𝐺𝑥(𝑘)+𝐻𝑣(𝑘),𝑥(𝑘+2)=𝐺𝑥(𝑘+1)+=𝐻𝑣(𝑘+1)𝐺²𝑥(𝑘)+𝐺𝐻𝑣(𝑘)+⋮𝐻𝑣(𝑘+1),𝑥(𝑘+𝑞)=𝐺𝑥(𝑘+𝑞−1)+=𝐻𝑣(𝑘+𝑞−1)𝐺^𝑞𝑥(𝑘)+𝐺^𝑞−1𝐻𝑣(𝑘)+𝐺^𝑞−2𝐻𝑣(𝑘+1)+⋯+𝐺𝐻𝑣(𝑘+𝑞−2)+𝐻𝑣(𝑘+𝑞−1). (27) Then, the output equation (1b) at time step 𝑘+𝑞 is given by𝑦(𝑘+𝑞)=𝐶𝑥(𝑘+𝑞)=𝐶𝐺^𝑞𝑥(𝑘)+𝐶𝐺^𝑞−1𝐻𝑣(𝑘)+⋯+𝐶𝐺𝐻𝑣(𝑘+𝑞−2)+𝐶𝐻𝑣(𝑘+𝑞−1). (28) Hence, these output equations, for 𝑘=0,1,…,𝑁−𝑞 with 𝑁>𝑞(𝑚+𝑝+1)−1, can be expressed as𝑦=𝐶𝐺^𝑞𝑥+𝜗𝑣, (29a) where𝑦=𝑦(𝑞)𝑦(𝑞+1)⋯𝑦(𝑁),,𝐶𝑥=𝑥(0)𝑥(1)⋯𝑥(𝑁−𝑞)𝜗=𝐻𝐶𝐺𝐻⋯𝐶𝐺^(𝑞−1)𝐻,𝑣=.𝑣(𝑞−1)𝑣(𝑞)⋯𝑣(𝑁−1)𝑣(𝑞−2)𝑣(𝑞−1)⋯𝑣(𝑁−2)⋮⋮⋱⋮𝑣(0)𝑣(1)⋯𝑣(𝑁−𝑞) (29b)

Note that the first term on the right-hand side of (29a) represents the effect of the preceding time steps, where 𝐺^𝑞 is sufficiently small and all the states in 𝑥 are bounded; so (29a) can be approximated by neglecting the first term on the right-hand side, such that𝑦≅𝜗𝑣. (30) It has the following least-squares solution:𝜗=𝑦𝑣⁺, (31) where 𝑣⁺ is the pseudoinverse of the matrix 𝑣 with 𝑣𝑣⁺=𝐼_𝑞. It clearly shows that once information of inputs and outputs is obtained, the observer Markov parameters 𝑌_𝑘=𝐶𝐺^𝑘−1𝐻, for 𝑘=1,2,3,…, can be known.

The Hankel matrix 𝐻_𝑎(𝑘−1) [3] is defined by𝐻_𝑎≜(𝑘−1)𝑌_𝑘𝑌_𝑘+1⋯𝑌_{𝑘+𝛽−1}𝑌_𝑘+1𝑌_𝑘+2⋯𝑌_𝑘+𝛽⋮⋮⋱⋮𝑌_{𝑘+𝛼−1}𝑌_𝑘+𝛼⋯𝑌_{𝑘+𝛼+𝛽−2},𝑘=1,2,3,…, (32) where 𝛼 and 𝛽 are sufficiently large positive integers. Furthermore, 𝐻_𝑎(𝑘−1) can be easily decomposed into𝐻_𝑎(𝑘−1)=Θ_𝐻𝑎𝐺^𝑘−1Φ_𝐻𝑎, (33) where Θ_𝐻𝑎 is the observability matrix defined byΘ_𝐻𝑎=𝐶^𝑇(𝐶𝐺)^𝑇(𝐶𝐺²)^𝑇⋯(𝐶𝐺^𝛼−1)^𝑇𝑇, (34a) and Φ_𝐻𝑎 is the controllability matrix given byΦ_𝐻𝑎=𝐻𝐺𝐻𝐺²𝐻⋯𝐺^𝛽−1𝐻. (34b)From (33), one knows that 𝐻_𝑎(0)=Θ_𝐻𝑎Φ_𝐻𝑎𝐻,for𝑘=1,_𝑎(1)=Θ_𝐻𝑎𝐺Φ_𝐻𝑎,for𝑘=2. (35a)(35b)When the observer Markov parameters are determined by using (31), the eigensystem realization algorithm (ERA) method is used to obtain the system realization estimates {𝐺,𝐻,𝐶} and the observer gain 𝐹 through the singular value decomposition (SVD) of the Hankel matrix. The ERA processes the SVD factorization of the block data matrix 𝐻_𝑎(𝑘−1), started for 𝑘=1, that is,𝐻_𝑎(0)=𝑉Σ𝑆^𝑇. (36) Here 𝑉 and 𝑆 are orthogonal matrices and Σ is a block diagonal matrix of the following form:ΣΣ=^𝑛000, (37) where Σ_𝑛=diag{𝜎₁𝜎₂⋯𝜎_𝑛min𝜎_𝑛min+1⋯𝜎_𝑛} contains monotonically nonincreasing entries 𝜎₁≥𝜎₂≥⋯≥𝜎_𝑛min≥𝜎_𝑛min+1≥⋯𝜎_𝑛≥0. Here some singular values 𝜎_𝑛min+1,…,𝜎_𝑛 are relatively small (𝜎_𝑛min+1≪𝜎_𝑛min) and negligible in the sense that they contain more noise information than system information. In order to construct the low-order observer of the system, defineΣ_𝑛min=𝜎₁0⋯000𝜎₂⋯00⋮⋮⋱⋮⋮00⋯𝜎_𝑛min−1000⋯0𝜎_𝑛min. (38) Hence, (36) can be rewritten as𝐻_𝑎(0)=𝑉_𝑛minΣ_𝑛min𝑆^𝑇_𝑛min=(𝑉_𝑛minΣ_𝑛^1/2_min)(Σ_𝑛^1/2_min𝑆^𝑇_𝑛min)=Θ_𝐻𝑎Φ_𝐻𝑎, (39) where 𝑉_𝑛min and 𝑆_𝑛min are composed of the first 𝑛_min columns of 𝑉 and 𝑆, respectively.

Similarly, applying (39) to (35b) gets𝐻_𝑎(1)=Θ_𝐻𝑎𝐺Φ_𝐻𝑎=𝑉_𝑛minΣ_𝑛^1/2_min𝐺Σ_𝑛^1/2_min𝑆^𝑇_𝑛min. (40) Hence, one obtains𝐺=Σ_𝑛^−1/2_min𝑉^𝑇_𝑛min𝐻_𝑎(1)𝑆_𝑛minΣ_𝑛^−1/2_min. (41) In other words, the reduced model of order 𝑛_min after deleting singular values 𝜎_𝑛min+1,…,𝜎_𝑛 is then considered as the robustly controllable and observable part of the realized system with an acceptable closed-loop performance. Observing (26b), (26c), (34a)–(34b), (39), and (41), the system realization estimates {𝐺,𝐻,𝐶} and the observer gain Fcan be found out as follows:𝐻=thefirst𝑚columnsofΣ_𝑛^1/2_min𝑆^𝑇_𝑛min,−𝐹=the(𝑚+1)thcolumntothe(𝑚+𝑝)thcolumnofΣ_𝑛^1/2_min𝑆^𝑇_𝑛min,𝐶=thefirst𝑝rowsof𝑉_𝑛minΣ_𝑛^1/2_min,𝐺=Σ_𝑛^−1/2_min𝑉^𝑇_𝑛min𝐻_𝑎(1)𝑆_𝑛minΣ_𝑛^−1/2_min−𝐹𝐶. (42) For system identification, SVD is very useful in determining the system order. In practice, the primary purpose of applying the OKID method is that the constructed observer satisfies the least-squares solution or acts the input-output map the same as a Kalman filter. If the data length is sufficiently long and the order of the observer is sufficiently large, the truncation error is negligible.

3.3. The Optimal Linearization Method

Many useful techniques for analysis and design of linear systems have been successfully developed and reported in literatures. Also, these techniques are easier to carry out performance analysis and controller design for linear systems than those for nonlinear systems. Unfortunately, most physical systems are nonlinear in nature. Therefore, it is desirable to find an effective linearization model for analyzing the dynamic characteristics of nonlinear systems [16]. The optimal linearization was first proposed in literature [17] for continuous-time nonlinear systems followed by stabilizing controller design for uncertain nonlinear systems using fuzzy models. The proposed optimal linearization at the operating state, not necessarily the equilibrium state, yields the exact linear model. Also it yields the optimal linear model defined by some convex constraint optimization criterion in the vicinity of the operating state. Taylor expansion is a common approach used; however for linearization, a truncated Taylor expansion usually results in an affine rather than linear model due to the generally nonvanishing constant term.

Consider the class of nonlinear systems governed by 𝑥(𝑘+1)=𝑓(𝑥(𝑘))+𝑞(𝑥(𝑘))𝑢(𝑘), (43) where 𝑓∶ℜ^𝑛→ℜ^𝑛 and 𝑞∶ℜ^𝑛→ℜ^𝑛×𝑚 are nonlinear with continuous partial derivatives with respect to the state variable at all steps 𝑘, where 𝑥(𝑘)∈ℜ^𝑛 is the state vector at time index 𝑘, and 𝑢(𝑘)∈ℜ^𝑚 is the control input vector at time index 𝑘. It is desired to have an exact local linearization system at an operating state of interest, say 𝑥_op(𝑘); namely,𝑥(𝑘+1)=𝐺(𝑘)𝑥(𝑘)+𝐻(𝑘)𝑢(𝑘), (44) where 𝐺(𝑘) and 𝐻(𝑘) are matrices of appropriate dimensions depending on 𝑥_op(𝑘). The linearization of the nonlinear system (43) is commonly represented by the truncated Taylor expansion as𝑥(𝑘+1)=𝑓(𝑥_eq(𝑘))+𝑞(𝑥_eq(𝑘))𝑢_eq(𝑘)+𝐴(𝑘)[𝑥(𝑘)−𝑥_eq(𝑘)]+𝐵(𝑘)[𝑢(𝑘)−𝑢_eq(𝑘)], (45) where {𝑥_eq(𝑘),𝑢_eq(𝑘)} is an equilibrium point, 𝐴(𝑘)=∇_𝑥𝑓(𝑥_eq(𝑘)), and 𝐵(𝑘)=∇_𝑥𝑞(𝑥_eq(𝑘)). ∇_𝑥𝑓(𝑥_eq(𝑘))∶ℜ^𝑛→ℜ^𝑛×𝑛 is the gradient of 𝑓(𝑥(𝑘)) (evaluated at 𝑥(𝑘)=𝑥_eq(𝑘)) with respect to 𝑥(𝑘). One can then represent (45) in the following form:𝑥(𝑘+1)=𝐴(𝑘)𝑥(𝑘)+𝐵(𝑘)𝑢(𝑘)−[𝐴(𝑘)𝑥_eq(𝑘)+𝐵(𝑘)𝑢_eq(𝑘)]. (46) Clearly, (46) is an affine rather than linear model due to the generally nonvanished term 𝐴(𝑘)𝑥_eq(𝑘)+𝐵(𝑘)𝑢_eq(𝑘). One exception is the trivial case where the equilibrium is zero, (𝑥_eq(𝑘),𝑢_eq(𝑘))=(0,0), which cannot, however, be ensured throughout a nonlinear control process.

Suppose that 𝑥_op(𝑘)≠0 is the given operating state; it is not necessarily an equilibrium of the given system (43). The goal is to construct a local linear model in terms of 𝑥(𝑘) and 𝑢(𝑘), so that it can well approximate the dynamical behaviors of the system in the vicinity of the operating state 𝑥_op(𝑘). In other words, one has𝑓(𝑥(𝑘))+𝑞(𝑥(𝑘))𝑢(𝑘)≈𝐺(𝑘)𝑥(𝑘)+𝐻(𝑘)𝑢(𝑘),forany𝑢(𝑘),𝑓(𝑥_op(𝑘))+𝑞(𝑥_op(𝑘))𝑢(𝑘)=𝐺(𝑘)𝑥_op(𝑘)+𝐻(𝑘)𝑢(𝑘),forany𝑢(𝑘). (47) Since the designed control input 𝑢(𝑘) is arbitrary, the following equality must be held:𝑞(𝑥_op(𝑘))=𝐻(𝑘), (48) so that (47) becomes quite simple; that is,𝑓(𝑥(𝑘))≈𝐺(𝑘)𝑥(𝑘),𝑓(𝑥_op(𝑘))=𝐺(𝑘)𝑥_op(𝑘). (49) To satisfy these, let 𝑔^𝑇_𝑖 denote the 𝑖th row of the matrix 𝐺(𝑘). Then (49) can be represented as𝑓_𝑖(𝑥(𝑘))≈𝑔^𝑇_𝑖𝑓𝑥(𝑘),𝑖=1,2,…,𝑛,_𝑖(𝑥_op(𝑘))=𝑔^𝑇_𝑖𝑥_op(𝑘),𝑖=1,2,…,𝑛, (50)(51) where 𝑓_𝑖∶ℜ^𝑛→ℜ is the 𝑖th component of 𝑓. Next, expanding the left-hand side of (50) about 𝑥_op(𝑘) and neglecting the second- and higher-order terms, one has𝑓_𝑖(𝑥_op(𝑘))+[∇_𝑥𝑓_𝑖(𝑥_op(𝑘))]^𝑇(𝑥(𝑘)−𝑥_op(𝑘))≈𝑔^𝑇_𝑖𝑥(𝑘). (52) Now, using (51), one can rewrite (52) as[∇_𝑥𝑓_𝑖(𝑥_op(𝑘))]^𝑇(𝑥(𝑘)−𝑥_op(𝑘))≈𝑔^𝑇_𝑖(𝑥(𝑘)−𝑥_op(𝑘)), (53) in which 𝑥(𝑘) is arbitrary but should be “close” to 𝑥_op(𝑘) so that the approximation is good. To determine a constant vector 𝑔^𝑇_𝑖, such that it is “as close as possible” to [∇_𝑥𝑓_𝑖(𝑥_op(𝑘))]^𝑇 and also satisfies 𝑔^𝑇_𝑖𝑥_op(𝑘)=𝑓_𝑖(𝑥_op(𝑘)), one may consider the following constrained minimization problem:min_𝑔𝑖1Ψ∶=2‖∇_𝑥𝑓_𝑖(𝑥_op(𝑘))−𝑔_𝑖‖²₂ (54) subject to𝑔^𝑇_𝑖𝑥_op(𝑘)=𝑓_𝑖(𝑥_op(𝑘)). (55) Notice that this is a convex constrained optimization problem; therefore, the first-order necessary condition for a minimum of Ψ is also sufficient, which is ∇_𝑔𝑖Ψ+𝜂∇_𝑔𝑖{𝑔^𝑇_𝑖𝑥_op(𝑘)−𝑓_𝑖(𝑥_op(𝑘))}=0, (56) where 𝜂 is the Lagrange multiplier. It follows from (56) that 𝑔_𝑖−∇_𝑥𝑓_𝑖(𝑥_op(𝑘))+𝜂𝑥_op(𝑘)=0. (57) Left multiplying (57) by 𝑥^𝑇_op(𝑘), one obtains𝑥𝜂=^𝑇op(𝑘)∇_𝑥𝑓_𝑖(𝑥_op(𝑘))−𝑓_𝑖(𝑥_op(𝑘))‖𝑥_op(𝑘)‖²₂,for𝑥_op(𝑘)≠0. (58) Substituting (58) into (57) gives𝑔_𝑖=∇_𝑥𝑓_𝑖(𝑥_op𝑓(𝑘))+^𝑖(𝑥_op(𝑘))−𝑥^𝑇_op(𝑘)∇_𝑥𝑓_𝑖(𝑥_op(𝑘))‖𝑥_op(𝑘)‖²₂𝑥_op(𝑘),for𝑥_op(𝑘)≠0. (59a) It is easily verified that when 𝑥_op(𝑘)=0, (57) yields 𝑔_𝑖=∇_𝑥𝑓_𝑖(0),for𝑥_op(𝑘)=0. (59b)

The controllability matrix Φ for the nonlinear system (43) at the operating state 𝑥_op(𝑘) is derived from the optimal linear model (𝐺(𝑘),𝐻(𝑘)) defined by (44), resulting inΦ=𝐻(𝑘)𝐺(𝑘)𝐻(𝑘)𝐺²(𝑘)𝐻(𝑘)⋯𝐺^𝑛−1(𝑘)𝐻(𝑘), (60) where 𝐺(𝑘) and 𝐻(𝑘) are constructed via the following rule: the 𝑖th columns of 𝐺(𝑘) and 𝐻(𝑘) are set to be zero whenever the 𝑖th components of 𝑥_op(𝑘) and of 𝑢(𝑘) are zero, respectively.

3.4. The Fast Online Identification and Observer of Unknown Nonlinear Stochastic Hybrid Systems

Through the optimal linearization methodology, the identified NARMAX model can be explicitly presented as the ARMAX model-based observer without any approximation at each sampling instant. Besides, based on OKID method, an initial parameter 𝜃_𝑖(0) of RELS algorithm (12) for STC will be presented in this section.

3.4.1. The Online Identification and Observer Based on NARMAX Model through Optimal Linearization

By the proposed NARMAX model shown in (11) for the STC, the discrete-time state-space innovation model (19) is constructed to design the control input 𝑢(𝑘) for simply controlling the unknown true nonlinear system. The elements 𝐺_𝑜, 𝐻_𝑜, and 𝐶_𝑜 of the control gains 𝐾_𝑑(𝑘) and 𝐸_𝑑(𝑘) could be obtained from the discrete-time state-space innovation model. Since the NARMAX model is nonlinear, the discrete-time state-space innovation model is linear; for this reason the optimal linearization method will be applied to the NARMAX model, so as to obtain a linear ARMAX model at operating states of interest without any approximation. Besides, it is also the optimal one in the sense of minimizing the optimization problem given in (54).

In this paper, a suitable NARMAX model with 𝑚-inputs, 𝑝-outputs and, 𝜌-time-steps could be selected as follows:𝑦_𝑖(𝑘)=𝐹_𝑖(𝑌)+𝐵_𝑖(𝑌)𝑈+𝐷_𝑖(𝑌)𝐸+𝜀_𝑖(𝑘),for𝑖=1,2,…,𝑝, (61) where𝑦(𝑘−𝑙)=𝑦₁(𝑘−𝑙)𝑦₂(𝑘−𝑙)⋯𝑦_𝑝(𝑘−𝑙)^𝑇𝑦,for𝑙=1,2,…,𝜌,𝑌=^𝑇(𝑘−1)𝑦^𝑇(𝑘−2)⋯𝑦^𝑇(𝑘−𝜌)^𝑇∈ℜ^𝜌𝑝×1𝐵_𝑖,(𝑌)=[𝑏_𝑖1(𝑌),…,𝑏_𝑖𝑚(𝑌),𝑏_𝑖(𝑚+1)(𝑌),…,𝑏_𝑖(2𝑚)(𝑌),…,𝑏_{𝑖,(𝜌−1)×𝑚+1}(𝑌),…,𝑏_{𝑖(𝜌×𝑚)}(𝑌)]∈ℜ^1×𝜌𝑚𝐷_𝑖,(𝑌)=[𝑑_𝑖1(𝑌),…,𝑑_𝑖𝑝(𝑌),𝑑_𝑖(𝑝+1)(𝑌),…,𝑑_𝑖(2𝑝)(𝑌),…,𝑑_{𝑖,(𝜌−1)×𝑝+1}(𝑌),…,𝑑_{𝑖(𝜌×𝑝)}(𝑌)]∈ℜ^1×𝜌𝑝𝑢,𝑢(𝑘−𝑙)=₁(𝑘−𝑙)𝑢₂(𝑘−𝑙)⋯𝑢_𝑚(𝑘−𝑙)^𝑇∈ℜ^𝑚×1𝑢,𝑈=^𝑇(𝑘−1)𝑢^𝑇(𝑘−2)⋯𝑢^𝑇(𝑘−𝜌)^𝑇∈ℜ^𝜌𝑚×1𝜀,𝜀(𝑘−𝑙)=₁(𝑘−𝑙)𝜀₂(𝑘−𝑙)⋯𝜀_𝑝(𝑘−𝑙)^𝑇∈ℜ^𝑝×1𝜀,𝐸=^𝑇(𝑘−1)𝜀^𝑇(𝑘−2)⋯𝜀^𝑇(𝑘−𝜌)^𝑇∈ℜ^𝜌𝑝×1. (62)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14)(15)(16)(17)(18) The prediction errors (𝜀_𝑖(𝑘−1),𝜀_𝑖(𝑘−2),…,𝜀_𝑖(0)) all have been identified from the NARMAX model based on the given control inputs (𝑢_𝑗(𝑘−1),𝑢_𝑗(𝑘−2),…,𝑢_𝑗(0)) and the output measurement data (𝑦_𝑖(𝑘),𝑦_𝑖(𝑘−1),…,𝑦_𝑖(0)) for 𝑖=1,2,…,𝑝 and 𝑗=1,2,…,𝑚. This chosen NARMAX model (61) has the advantage that 𝐵_𝑖(𝑌)𝑈 and 𝐷_𝑖(𝑌)𝐸 are linear in all items of 𝑢 and 𝜀, respectively.

To separate linear/nonlinear variables, this NARMAX model (61) can be rewritten as a linear-in-the-parameter expression as follows:𝑦_𝑖(𝑘)=𝐹_𝑖(𝑌)+𝐵_𝑖(𝑌)𝑈+𝐷_𝑖(𝑌)𝐸+𝜀_𝑖=(𝑘)_{𝜉𝑖𝑗=1}𝜃_𝑖𝑗𝜙_𝑖𝑗(𝑘)+𝜀_𝑖(𝑘)=𝜃^𝑇_𝑖𝜙_𝑖(𝑘)+𝜀_𝑖(𝑘), (63) where𝜃_𝑖=𝜃_𝑖1𝜃_𝑖2⋯𝜃_𝑖𝜉𝑖^𝑇,𝜙_𝑖(𝑘)=𝜙_𝑖1(𝑘)𝜙_𝑖2(𝑘)⋯𝜙_𝑖𝜉𝑖(𝑘)^𝑇. (64)

Estimate the parameter vector 𝜃_𝑖 by the RELS algorithm to make the output of the model approximate the true system output. Although 𝐹_𝑖(𝑌) only includes output delays, it is a nonlinear function of 𝑦. Therefore, 𝐹_𝑖(𝑌) needs to be linearized by the optimal linearization method in order to find the suboptimal controller.

Performing the optimal linearization approach on 𝐹_𝑖(𝑌) gives𝐹_𝑖(𝑌)=𝐴_𝑖𝑌, (65a) where 𝐴^𝑇_𝑖=(∇_𝑌𝐹_𝑖(𝐹𝑌)+^𝑖(𝑌)−𝑌^𝑇∇_𝑌𝐹_𝑖(𝑌)‖𝑌‖²₂𝑌)|_𝑌=𝑌_op≜𝑎_𝑖1⋯𝑎_𝑖𝑝𝑎_𝑖(𝑝+1)⋯𝑎_{𝑖(2×𝑝)}⋯𝑎_{𝑖,(𝜌−1)×𝑝+1}⋯𝑎_{𝑖(𝜌×𝑝)}^𝑇∈ℜ^𝜌𝑝×1. (65b)Substituting (65a) into (63) with𝐵_𝑖=𝐵_𝑖(𝑌)|_𝑌=𝑌_op≜[𝑏̃_𝑖1,…,𝑏̃_𝑖𝑚,𝑏̃_𝑖(𝑚+1),…,𝑏̃_𝑖(2𝑚)𝑏̃,…,_{𝑖,(𝜌−1)×𝑚+1},…,𝑏̃_{𝑖(𝜌×𝑚)}𝐷],_𝑖=𝐷_𝑖(𝑌)|_𝑌=𝑌_op≜[𝑑_𝑖1,…,𝑑_𝑖𝑝,𝑑_𝑖(𝑝+1),…,𝑑_𝑖(2𝑝)𝑑,…,_{𝑖,(𝜌−1)×𝑝+1},…,𝑑_{𝑖(𝜌×𝑝)}] (66)gives 𝑦_𝑖(𝑘)−𝐴_𝑖𝑌=𝐵_𝑖𝑈+𝐷_𝑖𝐸+𝜀_𝑖(𝑘). (67) Let 𝑧⁻¹ denote the backward shift operator. Equation (67) can be rewritten as𝑦_𝑖(𝑘)−𝑎_𝑖1𝑧⁻¹𝑦₁(𝑘)−⋯−𝑎_𝑖𝑝𝑧⁻¹𝑦_𝑝(𝑘)−𝑎_𝑖(𝑝+1)𝑧⁻²𝑦₁𝑎(𝑘)+⋯−_𝑖(2𝑝)𝑧⁻²𝑦_𝑝(𝑘)−⋯−𝑎_{𝑖,(𝜌−1)×𝑝+1}𝑧^−𝜌𝑦₁𝑎(𝑘)−⋯−_{𝑖(𝜌×𝑝)}𝑧^−𝜌𝑦_𝑝=𝑏̃(𝑘)_𝑖1𝑧⁻¹𝑢₁(𝑘)+⋯+𝑏̃_𝑖𝑚𝑧⁻¹𝑢_𝑚(𝑘)+𝑏̃_𝑖(𝑚+1)𝑧⁻²𝑢₁𝑏̃(𝑘)+⋯+_{𝑖(2×𝑚)}𝑧⁻²𝑢_𝑚(𝑘)+⋯+𝑏̃_{𝑖,(𝜌−1)×𝑚+1}𝑧^−𝜌𝑢₁𝑏̃(𝑘)+⋯+_{𝑖(𝜌×𝑚)}𝑧^−𝜌𝑢_𝑚(𝑘)+𝜀_𝑖(𝑘)+𝑑_𝑖1𝑧⁻¹𝜀₁𝑑(𝑘)+⋯+_𝑖𝑝𝑧⁻¹𝜀_𝑝(𝑘)+𝑑_𝑖(𝑝+1)𝑧⁻²𝜀₁𝑑(𝑘)+⋯+_{𝑖(2×𝑝)}𝑧⁻²𝜀_𝑝(𝑘)+⋯+𝑑_{𝑖,(𝜌−1)×𝑝+1}𝑧^−𝜌𝜀₁𝑑(𝑘)+⋯+_{𝑖(𝜌×𝑝)}𝑧^−𝜌𝜀_𝑝(𝑘),for𝑖=1,2,…,𝑝. (68)

After combining (68) for 𝑖=1,2,…,𝑝, a disturbed output of the plant in the ARMAX model takes the following form:𝐺_𝐿(𝑧⁻¹)𝑦(𝑘)=𝐻_𝐿(𝑧⁻¹)𝑢(𝑘)+𝐷_𝐿(𝑧⁻¹)𝜀(𝑘), (69a) where 𝐺_𝐿(𝑧⁻¹)=𝐼_𝑝+𝐺_𝑜1𝑧⁻¹+⋯+𝐺_𝑜𝜌𝑧^−𝜌𝐺_𝑜𝑙,(𝑙=1,2,…,𝜌)∈ℜ^𝑝×𝑝𝐻_𝐿,(𝑧⁻¹)=𝐻_𝑜1𝑧⁻¹+𝐻_𝑜2𝑧⁻²+⋯+𝐻_𝑜𝜌𝑧^−𝜌𝐻_𝑜𝑙,(𝑙=1,2,…,𝜌)∈ℜ^𝑝×𝑚𝐷_𝐿,(𝑧⁻¹)=𝐼_𝑝+𝐷_𝑜1𝑧⁻¹+⋯+𝐷_𝑜𝜌𝑧^−𝜌𝐷_𝑜𝑙,(𝑙=1,2,…,𝜌)∈ℜ^𝑝×𝑝𝐺_𝑜𝑙,=−𝑎_{1,(𝑙−1)×𝑝+1}−𝑎_{1,(𝑙−1)×𝑝+2}⋯−𝑎_{1(𝑙×𝑝)}−𝑎_{2,(𝑙−1)×𝑝+1}−𝑎_{2,(𝑙−1)×𝑝+2}⋯−𝑎_{2(𝑙×𝑝)}−𝑎⋮⋮⋱⋮_{𝑝,(𝑙−1)×𝑝+1}−𝑎_{𝑝,(𝑙−1)×𝑝+2}⋯−𝑎_{𝑝(𝑙×𝑝)}𝐻_𝑜𝑙,=𝑏̃_{1,(𝑙−1)×𝑚+1}𝑏̃_{1,(𝑙−1)×𝑚+2}⋯𝑏̃_{1(𝑙×𝑚)}𝑏̃_{2,(𝑙−1)×𝑚+1}𝑏̃_{2,(𝑙−1)×𝑚+2}⋯𝑏̃_{2(𝑙×𝑚)}𝑏̃⋮⋮⋱⋮_{𝑝,(𝑙−1)×𝑚+1}𝑏̃_{𝑝,(𝑙−1)×𝑚+2}⋯𝑏̃_{𝑝(𝑙×𝑚)}𝐷_𝑜𝑙,=𝑑_{1,(𝑙−1)×𝑝+1}𝑑_{1,(𝑙−1)×𝑝+2}⋯𝑑_{1(𝑙×𝑝)}𝑑_{2,(𝑙−1)×𝑝+1}𝑑_{2,(𝑙−1)×𝑝+2}⋯𝑑_{2(𝑙×𝑝)}𝑑⋮⋮⋱⋮_{𝑝,(𝑙−1)×𝑝+1}𝑑_{𝑝,(𝑙−1)×𝑝+2}⋯𝑑_{𝑝(𝑙×𝑝)}. (69b)(19)(20)(21)(22)(23)(24)(25)(26)

An alternative representation of the ARMAX model (69a) and (69b) is given by𝑦(𝑘)=𝐺_𝐿⁻¹(𝑧⁻¹)𝐻_𝐿(𝑧⁻¹)𝑢(𝑘)+𝐺_𝐿⁻¹(𝑧⁻¹)𝐷_𝐿(𝑧⁻¹)𝜀(𝑘), (70) which is in the left matrix fraction description form (LMFD) [6]. The first and second terms of (70) share the same left characteristic matrix polynomial 𝐺_𝐿⁻¹(𝑧⁻¹), which represents the effects of control inputs and noises. Once 𝐺_𝐿⁻¹(𝑧⁻¹) has been specified to characterize the dynamics of the plant, the error vector model 𝐺_𝐿⁻¹(𝑧⁻¹)𝐷_𝐿(𝑧⁻¹)𝜀(𝑘) presents an adjustable moving average process of 𝐷_𝐿(𝑧⁻¹)𝜀(𝑘) regarded as the colored noise input. Under the linearized ARMAX model, a system in an observable block companion form can be described by the following state-space innovation form [18–20]:𝑥_𝑜(𝑘+1)=𝐺_𝑜𝑥_𝑜(𝑘)+𝐻_𝑜𝑢(𝑘)+𝐾_𝑜(𝑘)𝑒_𝑜𝑒(𝑘),_𝑜(𝑘)=𝑦(𝑘)−𝐶_𝑜𝑥_𝑜(𝑘), (71a)(71b) where𝐺_𝑜=−𝐺_𝑜1𝐼_𝑝0_𝑝⋯0_𝑝−𝐺_𝑜20_𝑝𝐼_𝑝⋯0_𝑝⋮⋮⋮⋱⋮−𝐺_{𝑜(𝜌−1)}0_𝑝0_𝑝⋯𝐼_𝑝−𝐺_𝑜𝜌0_𝑝0_𝑝⋯0_𝑝𝐻_𝑜,=𝐻^𝑇_𝑜1𝐻^𝑇_𝑜2⋯𝐻^𝑇_𝑜𝜌^𝑇𝐶_𝑜,=𝐼_𝑝0_𝑝⋯0_𝑝𝐾_𝑜,=𝐾_𝑜1𝐾_𝑜2⋯𝐾_𝑜𝜌,𝐾_𝑜𝑙=𝐷_𝑜𝑙−𝐺_𝑜𝑙𝑥for𝑙=1,2,…,𝜌,_𝑜(𝑘)=𝑥^𝑇_𝑜1(𝑘)𝑥^𝑇_𝑜2(𝑘)⋯𝑥^𝑇_𝑜𝜌(𝑘)^𝑇,𝑥_𝑜𝑙(𝑘)∈ℜ^𝑝for𝑙=1,2,…,𝜌. (71c)(71d)(71e)(71f)(71g)𝑥_𝑜(𝑘) is the estimate of the system state 𝑥_𝑜(𝑘) in the observer coordinates. In this process, the matrices 𝐺_𝑜𝑙, 𝐻_𝑜𝑙, and 𝐷_𝑜𝑙 are calculated as 𝐺_𝑜𝑙(𝑘), 𝐻_𝑜𝑙(𝑘), and 𝐷_𝑜𝑙(𝑘) at each sampling time 𝑘𝑇, respectively.

However, the zeros of 𝐷_𝐿(𝑧⁻¹) in (69a) and (69b) may not be all in the unit circle centered at the origin, so that the eigenvalues of the Kalman filter (71a)–(71g) with gain 𝐾_𝑜(𝑘) may not all lie in the unit circle centered at the origin, either. Instead of the Kalman gain 𝐾_𝑜(𝑘), one could design the digital estimator gain 𝐿_𝑜(𝑘) to replace 𝐾_𝑜(𝑘). The digital estimator gain 𝐿_𝑜(𝑘) will be indirectly designed by employing the discrete-time observer design based on the estimated process model, rather than indirectly estimated from the identified parameter of the NARMAX model. If (𝐺^𝑇_𝑜(𝑘),𝐶^𝑇_𝑜) is stabilizable, there exists a digital estimator gain 𝐿_𝑜(𝑘) so that the closed-loop estimator matrix 𝐺_𝑜(𝑘)−𝐿_𝑜(𝑘)𝐶_𝑜 is asymptotically stable. The estimator gain is𝐿^𝑇_𝑜(𝑘)=(𝐶_𝑜𝑃_𝑜(𝑘)𝐶^𝑇_𝑜+𝑅_𝑜)⁻¹𝐶_𝑜𝑃_𝑜(𝑘)𝐺^𝑇_𝑜(𝑘). (72a) Here 𝑃_𝑜(𝑘) is the solution to the following Riccati equation:𝐺_𝑜(𝑘)𝑃_𝑜(𝑘)𝐺^𝑇_𝑜(𝑘)−𝑃_𝑜(𝑘)−(𝐺_𝑜(𝑘)𝑃_𝑜(𝑘)𝐶^𝑇_𝑜)×(𝐶_𝑜𝑃_𝑜(𝑘)𝐶^𝑇_𝑜+𝑅_𝑜)⁻¹(𝐶_𝑜𝑃_𝑜(𝑘)𝐺^𝑇_𝑜(𝑘))+𝑄_𝑜=0, (72b)in which 𝑄_𝑜≥0 and 𝑅_𝑜>0 with appropriate dimensions. Then the state-space innovation form (71a) becomes𝑥_𝑜(𝑘+1)=𝐺_𝑜(𝑘)𝑥_𝑜(𝑘)+𝐻_𝑜(𝑘)𝑢(𝑘)+𝐿_𝑜(𝑘)𝑒_𝑜(𝑘). (73) Thus, one can design the digital controller𝑢(𝑘)=−𝐾_𝑑(𝑘)𝑥_𝑜(𝑘)+𝐸_𝑑(𝑘)Γ^∗(𝑘) (74) by the observer (73).

3.4.2. The Initial Parameters of NARMAX Model Based on OKID

The initial parameter vector 𝜃_𝑖(0) of NARMAX model will affect the convergent speed of RELS process. In order to increase the convergent speed of RELS algorithm, 𝜃_𝑖(0) will be figured out via the following method. The first step is to perform the offline system identification scheme by using the OKID method to obtain the discrete system realization estimates {𝐺,𝐻,𝐶} and observer gain 𝐹. Then, transfer these matrices into the corresponding observer form {𝐺_𝑜,𝐻_𝑜,𝐶_𝑜} and observer gain 𝐾_𝑜 (71a)–(71g) byΘ=(𝐶𝐺^𝜌−1)^𝑇(𝐶𝐺^𝜌−2)^𝑇⋯(𝐶𝐺)^𝑇𝐶^𝑇𝑇,𝑇_𝑜1=Θ⁻¹𝐼_𝑝0_𝑝⋯0_𝑝^𝑇∈ℜ^𝑛×𝑝,𝑇_𝑜=𝐺^𝜌−1𝑇_𝑜1𝐺^𝜌−2𝑇_𝑜1⋯𝐺𝑇_𝑜1𝑇_𝑜1∈ℜ^𝑛×𝑛,𝐺_𝑜=𝑇_𝑜⁻¹𝐺𝑇_𝑜∈ℜ^𝑛×𝑛,𝐻_𝑜=𝑇_𝑜⁻¹𝐻=𝐻^𝑇_𝑜1𝐻^𝑇_𝑜2⋯𝐻^𝑇_𝑜𝜌^𝑇∈ℜ^𝑛×𝑚,𝐶_𝑜=𝐶𝑇_𝑜=𝐼_𝑝0_𝑝⋯0_𝑝∈ℜ^𝑝×𝑛,𝐾_𝑜=𝑇_𝑜⁻¹𝐹=𝐾_𝑜1𝐾_𝑜2⋯𝐾_𝑜𝜌. (75) Based on (69a)-(69b) and (71a)–(71g) one could get an ARMAX model as follows:𝑦(𝑘)+𝐺_𝑜1𝑦(𝑘−1)+𝐺_𝑜2𝑦(𝑘−2)+⋯+𝐺_𝑜𝜌=𝐻𝑦(𝑘−𝜌)_𝑜1𝑢(𝑘−1)+𝐻_𝑜2𝑢(𝑘−2)+⋯+𝐻_𝑜𝜌+𝐷𝑢(𝑘−𝜌)_𝑜1𝜀(𝑘−1)+𝐷_𝑜2𝜀(𝑘−2)+⋯+𝐷_𝑜𝜌𝜀(𝑘−𝜌), (76) where 𝐷_𝑜𝑙=𝐺_𝑜𝑙+𝐾_𝑜𝑙,𝑙=1,2,…,𝜌.

Comparing (76) with (69a) and (69b) obtains the coefficient matrices {𝐴_𝑖,𝐵_𝑖,𝐷_𝑖} of the linearized ARMAX model defined by (65a)-(65b), (66), and (67). Rearranging the right-hand sides of (65b) and (66), respectively, gives𝐴_𝑖=[∇_𝑌𝐹_𝑖(𝑌_op𝐹)+^𝑖(𝑌_op)−𝑌^𝑇_op∇_𝑌𝐹_𝑖(𝑌_op)‖𝑌_op‖²₂𝑌_op]^𝑇=𝑎_𝑖1𝑎_𝑖2⋯𝑎_{𝑖𝑛𝐹𝑖}𝐹_𝑖(𝑌_op𝐵),_𝑖=𝐵_𝑖(𝑌)|_𝑌=𝑌_op=𝑏_𝑖1𝑏_𝑖2⋯𝑏_{𝑖𝑛𝐵𝑖}𝐵_𝑖(𝑌_op𝐷),_𝑖=𝐷_𝑖(𝑌)|_𝑌=𝑌_op=𝑑_𝑖1𝑑_𝑖2⋯𝑑_{𝑖𝑛𝐷𝑖}𝐷_𝑖(𝑌_op), (77) where 𝑖=1,2,…,𝑝. The unknown parameters 𝑎_𝑖𝑙, 𝑏_𝑖𝑚, and 𝑑_𝑖𝑛 are the undetermined coefficients of 𝐹_𝑖(𝑌),𝐵_𝑖(𝑌), and 𝐷_𝑖(𝑌) for 𝑙=1,2,…,𝑛_𝐹𝑖, 𝑚=1,2,…,𝑛_𝐵𝑖, and 𝑛=1,2,…,𝑛_𝐷𝑖, respectively. 𝐹_𝑖(𝑌_op), 𝐵_𝑖(𝑌_op), and 𝐷_𝑖(𝑌_op) are composed of the corresponding terms with respect to 𝑎_𝑖𝑙, 𝑏_𝑖𝑚, and 𝑑_𝑖𝑛, respectively.

Note that 𝜉_𝑖=𝑛_𝐹𝑖+𝑛_𝐵𝑖+𝑛_𝐷𝑖. These equations show the relationship between the unknown parameters of the NARMAX model (63) and the known parameters of the linearized ARMAX model (67). These unknown parameters can be solved by applying the pseudo inverse operation to (77). The solutions can be obtained as follows:𝑎_𝑖1𝑎_𝑖2⋯𝑎_{𝑖𝑛𝐹𝑖}=𝐴_𝑖𝐹⁺_𝑖(𝑌_op𝑏̂),_𝑖1𝑏̂_𝑖2⋯𝑏̂_{𝑖𝑛𝐵𝑖}=𝐵_𝑖𝐵⁺_𝑖(𝑌_op𝑑),_𝑖1𝑑_𝑖2⋯𝑑_{𝑖𝑛𝐷𝑖}=𝐷_𝑖𝐷⁺_𝑖(𝑌_op). (78) Consequently, the parameter vector estimate𝜃_𝑖=𝑎_𝑖1𝑎_𝑖2⋯𝑎_{𝑖𝑛𝐹𝑖}𝑏̂_𝑖1𝑏̂_𝑖2⋯𝑏̂_{𝑖𝑛𝐵𝑖}𝑑_𝑖1𝑑_𝑖2⋯𝑑_{𝑖𝑛𝐷𝑖}^𝑇=𝐴_𝑖𝐹⁺_𝑖(𝑌_op)𝐵_𝑖𝐵⁺_𝑖(𝑌_op)𝐷_𝑖𝐷⁺_𝑖(𝑌_op)^𝑇 (79) can be served as the initial parameter 𝜃_𝑖(0) of RELS algorithm (12) for STC.

4. NARMAX Model-Based State-Space Self-Tuning Control for Unknown Nonlinear Stochastic Hybrid Systems with Initial OKID-Estimated Parameters

Consider the class of continuous-time nonlinear stochastic systems as follows.𝑥̇(𝑡)=𝑓(𝑥(𝑡))+𝑞(𝑥(𝑡))𝑢(𝑡)+𝑤(𝑡),𝑦(𝑡)=ℎ(𝑥(𝑡))+𝑣(𝑡), (80) where 𝑓∶ℜ^𝑛→ℜ^𝑛, 𝑞∶ℜ^𝑛→ℜ^𝑛×𝑚, ℎ∶ℜ^𝑛→ℜ^𝑝, 𝑢(𝑘)∈ℜ^𝑚 is the control input, 𝑥(𝑘)∈ℜ^𝑛 is the state vector, 𝑦(𝑡)∈ℜ^𝑝 is the measurable output vector, and 𝑤(𝑡)∈ℜ^𝑛 and 𝑣(𝑡)∈ℜ^𝑝 are uncorrelated white noise processes. It is assumed that 𝑓(𝑥(𝑡)), 𝑞(𝑥(𝑡)), ℎ(𝑥(𝑡)), 𝑤(𝑡), and 𝑣(𝑡) are unknown, and 𝑥(𝑘) is inaccessible. The structure of the proposed state-space self-tuning control with the NARMAX model is shown in Figure 3.

Figure 3 

Structure of the hybrid state-space self-tuning control with NARMAX model and OKID for unknown nonlinear stochastic hybrid systems.

The design procedure is given as follows.

Step 1. For the unknown continuous-time nonlinear stochastic system (7a) and (7b), choose an appropriate NARMAX model (61) to be identified for self-tuning control.(i)Perform the offline system identification scheme described in Section 3.2 to obtain system and observer-gain Markov parameters of the OKID model; next, use the ERA method to obtain the system realization estimates 𝐺, 𝐻, 𝐶, and observer gain 𝐹 and then transfer them into observer form 𝐺_𝑜, 𝐻_𝑜, 𝐶_𝑜, and 𝐾_𝑜.(ii)based on the state-space innovation form (71a)–(71g) and the optimal linearization of NARMAX model in (65a)–(69b), the initial parameter 𝜃_𝑖(0) can be obtained by using (79).

Step 2. Identify the parameter of the NARMAX model at each sampling period T.(i)Let the number of 𝜃_𝑖 be 𝜉_𝑖. Also, set 𝑆_𝑖(0)=𝛾_𝑖𝐼_𝜉𝑖>0, 0<𝜆₀<1, 0.9<𝜆(0)<1, 𝑥_𝑜(0), 𝑘=1, and 𝜃_𝑖(0) obtained by OKID in Step 1.(ii)For online identifying the unknown continuous-time nonlinear stochastic system (80) with the piecewise constant control input, one utilizes the information of inputs and outputs to determine the updated parameter estimate 𝜃_𝑖(𝑘) at each sampling interval 𝑇 by RELS algorithm.

Step 3. Linearize the NARMAX model at each sampling period T:
Based on the updated parameter estimate 𝜃_𝑖(𝑘), linearizing the NARMAX model by the optimal linearization methodology and using (69a)-(69b) and (71a)–(71g) yields 𝐺_𝑜(𝑘), 𝐻_𝑜(𝑘), and 𝐶_𝑜.

Step 4. Estimate the predicted state at each sampling period T:
Select appropriate {𝑄_𝑜,𝑅_𝑜} in (72b) to have the high-gain property digital estimator gain:𝐿^𝑇_𝑜(𝑘)=(𝐶_𝑜𝑃_𝑜(𝑘)𝐶^𝑇_𝑜+𝑅_𝑜)⁻¹𝐶_𝑜𝑃_𝑜(𝑘)𝐺^𝑇_𝑜(𝑘), (81) where 𝑃_𝑜(𝑘) is the positive definite and symmetric solution to the Riccati equation (72b). The associated state-space innovation form is given by 𝑥_𝑜(𝑘∣𝑘−1)=𝐺_𝑜(𝑘−1)𝑥_𝑜(𝑘−1∣𝑘−2)+𝐻_𝑜(𝑘−1)𝑢(𝑘−1)+𝐿_𝑜(𝑘−1)𝑒_𝑜𝑒(𝑘−1),_𝑜(𝑘−1)=𝑦(𝑘−1)−𝐶_𝑜𝑥_𝑜(𝑘−1∣𝑘−2). (82)
From the above three equations, one can predict the state 𝑥_𝑜(𝑘∣𝑘−1).

Step 5. Generate the digital control input at each sampling period 𝑇.(i)Select appropriate weighting matrices {𝑄_𝑑,𝑅_𝑑} for the Riccati equation (5) such that the desired digital control law possesses the high-gain property.(ii)Compute the digital control gains 𝐾_𝑑(𝑘) and 𝐸_𝑑(𝑘) by the digital control formula (4) as follows: 𝐾_𝑑(𝑘)=[𝐻^𝑇_𝑜(𝑘)𝑃(𝑘)𝐻_𝑜(𝑘)+𝑅_𝑑]⁻¹𝐻^𝑇_𝑜(𝑘)𝑃(𝑘)𝐺_𝑜𝐸(𝑘),_𝑑(𝑘)=[𝐻^𝑇_𝑜(𝑘)𝑃(𝑘)𝐻_𝑜(𝑘)+𝑅_𝑑]⁻¹×𝐻^𝑇_𝑜(𝑘)[𝐼−(𝐺_𝑜(𝑘)−𝐻_𝑜(𝑘)𝐾_𝑑)^𝑇]⁻¹𝐶^𝑇_𝑜𝑄_𝑑, (83) where 𝑃(𝑘) is the positive definite symmetric solution of the following Riccati equation: 𝐺^𝑇_𝑜(𝑘)𝑃(𝑘)𝐺_𝑜(𝑘)−𝑃(𝑘)+𝐶^𝑇_𝑜𝑄_𝑜𝐶_𝑜−𝐺^𝑇_𝑜(𝑘)𝑃(𝑘)𝐻_𝑜(𝑘)×[𝐻^𝑇_𝑜(𝑘)𝑃(𝑘)𝐻_𝑜(𝑘)+𝑅_𝑜]⁻¹𝐻^𝑇_𝑜(𝑘)𝑃(𝑘)𝐺_𝑜(𝑘)=0. (84)(iii)Generate the digital PAM control input at each sampling period 𝑇 by using (73) and (74) with 𝐾_𝑑(𝑘) and 𝐸_𝑑(𝑘) obtained in Step 5(ii) 𝑢(𝑡)=𝑢(𝑘𝑇)=−𝐾_𝑑(𝑘𝑇)𝑥_𝑜(𝑘𝑇∣𝑘𝑇−𝑇)+𝐸_𝑑(𝑘𝑇)Γ(𝑘𝑇+𝑇),for𝑘𝑇≤𝑡<𝑘𝑇+𝑇. (85) Note that the reference input Γ(𝑡)=Γ(𝑘𝑇+𝑇), for 𝑘𝑇≤𝑡<𝑘𝑇+𝑇, is based on the tracking purpose.(iv)Calculate the PWM control input 𝑢_𝑤(𝑡) by using (9) and (10).(v)Set 𝑘=𝑘+1. Go to Step 2(ii) and continue the adaptive control process.

5. Self-Tuning Control with Fault Tolerance

5.1. Problem Statement

Consider again the nonlinear continuous-time stochastic system given by (80). If the system states and inputs are in partial faults, let the faulty system be characterized by𝑥̇(𝑡)=𝑓(𝑥)+𝑞(𝑥)𝑢(𝑡)+^𝜍_𝑣=1𝛽_𝑣(𝑡−𝜏_𝑣)[𝑓_𝑣(𝑥)+𝑞_𝑣(𝑢)]+𝑤(𝑡), (86) where 𝑓_𝑣(𝑥)+𝑞_𝑣(𝑢) represents the dynamic changes caused by the unknown and unanticipated failure modes of states and inputs. Two typical faults, gradual fault and abrupt fault, are considered online. Their characteristics are described by the time-varying function 𝛽_𝑣(𝑡) as shown in Figure 4 [21], where 𝑢_𝑠(𝑡) denotes the unit-step function of time 𝑡. 𝑓_𝑣(𝑥)+𝑞_𝑣(𝑢), 𝛽_𝑣(𝑡), and 𝜏_𝑣 are unknown due to the possible occurrence of unanticipated faults. If 𝜍=1, then a single-fault happened; 𝜍=2,3,4,… means the multiple-fault case.

(a)

(b)

Figure 4 

(a) The gradual failure function. (b) The abrupt failure function.

By the optimal linearization method described in Section 3.3, the nonlinear faulty system (86) can be accurately linearized as the following linear time-varying model:𝑥̇(𝑡)=[𝐺(𝑡)+^𝜍_𝑣=1𝛽_𝑣(𝑡−𝜏_𝑣)𝐺_𝑣(𝑡)]𝑥(𝑡)+[𝐻(𝑡)+^𝜍_𝑣=1𝛽_𝑣(𝑡−𝜏_𝑣)𝐻_𝑣(𝑡)]𝑢(𝑡)+𝑤(𝑡). (87) It is clear that the system could contain large uncertainties when the failure dynamics are large. Under this situation, the controller has to take an appropriate control action for the uncertainties occurring at any time instant 𝜏_𝑣. This is an adaptive control problem, in which controller parameters are adjusted according to the system parameter estimates. In this paper, the method based on the modified STC scheme is proposed to accomplish the fault-tolerant control.

Three assumptions required for the proposed fault-tolerant control are addressed as follows.

Assumption. The system is controllable and observable even if faults occur.

Assumption. The system is persistently excited by control inputs.

Assumption. Before the fault occurs, the system is healthy or has fully recovered from the previous fault.

5.2. Modified Self-Tuning Control with Active Fault Tolerance

The STC scheme should be modified to cope with parameter variations due to system faults. When a partial fault occurs, the system parameters vary accordingly. The time-varying parameter estimates obtained via the RELS algorithm in the conventional STC scheme would give large parameter errors and lead to poor system performances. However, based on the Kalman filter interpretation algorithm of RELS method [22], a modified scheme is proposed to estimate parameter variations. The modified state-space STC scheme can be applied to multivariable stochastic faulty systems without requiring prior knowledge of system parameters and noise properties.

In short, in the beginning, a healthy but unknown system is well tuned by the conventional STC scheme, and then the self-tuning structure with the reset covariance matrix of the parameter estimate is modified to enhance the parameter estimation and output response when the system states and/or inputs are partially faulty [9].

It postulates that the parameter vector is time-varying; that is, it is modeled as a random walk or a drift:𝜃_𝑖(𝑘+1)=𝜃_𝑖(𝑘)+𝑤_𝑖𝑦(𝑘),_𝑖(𝑘)=𝜃^𝑇_𝑖(𝑘)𝜙_𝑖(𝑘)+𝑣_𝑖(𝑘),𝐸(𝑤_𝑖(𝑘)𝑤^𝑇_𝑖(𝑘))=𝑄_𝑖,𝐸(𝑣_𝑖(𝑘)𝑣^𝑇_𝑖(𝑘))=𝑅_𝑖, (88a)(88b)(88c)(88d)where 𝑖=1,2,…,𝑝. {𝑤_𝑖(𝑘)} and {𝑣_𝑖(𝑘)} are white Gaussian noise sequences with zero means. Applying the Kalman filter to this model (88a) and (88d) gives the following recursive algorithm:𝜃_𝑖(𝑘)=𝜃_𝑖(𝑘−1)+𝑀_𝑖(𝑘)𝜀_𝑖𝑀(𝑘),_𝑖(𝑘)=𝑆_𝑖(𝑘−1)𝜙_𝑖(𝑘)𝑅_𝑖+𝜙^𝑇_𝑖(𝑘)𝑆_𝑖(𝑘−1)𝜙_𝑖(𝑘),𝑆_𝑖(𝑘)=𝑆_𝑖(𝑘−1)−𝑆_𝑖(𝑘−1)𝜙_𝑖(𝑘)𝜙^𝑇_𝑖(𝑘)𝑆_𝑖(𝑘−1)𝑅_𝑖+𝜙^𝑇_𝑖(𝑘)𝑆_𝑖(𝑘−1)𝜙_𝑖(𝑘)+𝑄_𝑖𝜀_𝑖,(𝑘)=𝑦_𝑖(𝑘)−𝜃^𝑇_𝑖(𝑘−1)𝜙_𝑖(𝑘), (89a)(89b)(89c)(89d) with𝑆_𝑖(0)=𝐸{[𝜃_𝑖(0)−𝜃_𝑖(0)][𝜃_𝑖(0)−𝜃_𝑖(0)]^𝑇}, (89e)where 𝑆_𝑖(𝑘) is the covariance matrix of the parameter estimate error. {𝜀_𝑖(𝑘)} will usually be (near) a white noise sequence corrupting the system output data if the model with the parameter estimate is in good agreement with its true system.

To modify Kalman filter interpretation of RELS method, some appropriate initializations of 𝑄_𝑖, 𝑅_𝑖, and 𝑆_𝑖(0) are assumed to be prespecified before the parameter estimation process. However, when unanticipated system faults occur during the process with large parameter variations, naturally they need to be reasonably reset.

To modify conventional STC process with the RELS estimation algorithm for the faulty system, we propose to approximate the initializations of 𝑄_𝑖, 𝑅_𝑖, and 𝑆_𝑖(0) to the following moving window-based statistical quantities: 𝑅_𝑖≈1𝑁−𝑘₁+1^𝑁_𝑘=𝑘1𝜀_𝑖(𝑘)𝜀^𝑇_𝑖(𝑘),𝑆_𝑖1(0)≈diag{diagonalelementsof𝑁−𝑘₁×+1_{𝑁𝑘=𝑘1}[(𝜃_𝑖(𝑘)−𝜃_𝑖(𝑁))(𝜃_𝑖(𝑘)−𝜃_𝑖(𝑁))^𝑇]}, (90)(91) where 𝑘₁ is the time index after the estimate 𝜃_𝑖(𝑘) is in steady state, and 𝑁 is the last time instant of the conventional STC. It should be noted that the elements of 𝜃_𝑖(𝑘) in (91) would not be independent with each other, when 𝑁−𝑘₁+1 is not large enough. Similarly, 𝑄_𝑖 can be approximated by comparing (88a) with (89a) as follows:𝑄_𝑖≈[𝑀_𝑖(𝑁)𝜀_𝑖(𝑁)][𝑀_𝑖(𝑁)𝜀_𝑖(𝑁)]^𝑇, (92) where𝑀_𝑖(𝑁)=𝑆_𝑖(0)𝜙_𝑖(𝑁)𝑅_𝑖+𝜙^𝑇_𝑖(𝑁)𝑆_𝑖(0)𝜙_𝑖(𝑁). (93)

The proposed STC with the algorithm in (89a)–(89e) and the initializations in (90)–(93) works only for the plant with slowly time-varying parameters. This can be interpreted by the fact that the initialized 𝑄_𝑖 obtained from (92), while the system is healthy, is so small in general; hence 𝜃_𝑖(𝑘)≅𝜃_𝑖(𝑘−1) in (89a). As a result, it cannot reflect the real parameter variations induced by the unanticipated system faults. Therefore, 𝑆_𝑖(𝑘), 𝑅_𝑖, and 𝑄_𝑖 in (89c) need to be appropriately reset when a fault is detected at time instant 𝑘_𝑓𝑇. Although the conventional STC with an appropriately reset forgetting factor 𝜆(𝑘_𝑓) could improve the estimation of parameter variations, the reset forgetting factor 𝜆(𝑘_𝑓) would need some trials for various failure modes. Nevertheless, the resets of 𝑆_𝑖(𝑘_𝑓), 𝑄_𝑖, and 𝑅_𝑖 proposed in [9] are a systematic approach for various failure modes.

Due to the fact that the parameter variations induced by faults are unknown, the rule of thumb to reset the covariance matrices of the parameter estimate 𝑆_𝑖(𝑘) is given as follows. The variation of the parameter estimate before and after the fault can be approximated as𝛿𝜃_𝑖(𝑘_𝑓)≈𝜃_𝑖(𝑘_𝑓)−𝜃_𝑖(𝑘_ℎ),for𝑘_ℎ<𝑘_𝑓, (94) where 𝜃_𝑖(𝑘_ℎ) is the parameter estimate of the healthy system and 𝜃_𝑖(𝑘_𝑓) is the parameter estimate of the faulty system. Then, based on the physical interpretation of (89e), 𝑆_𝑖(𝑘−1) in (89c) can be reasonably reset as𝑆_𝑖(𝑘_𝑓𝜃−1)≡diag{diagonalelementsof[_𝑖(𝑘_𝑓)−𝜃_𝑖(𝑘_ℎ𝜃)]×[_𝑖(𝑘_𝑓)−𝜃_𝑖(𝑘_ℎ)]^𝑇}≈𝛿²diag{diagonalelementsof𝜃_𝑖(𝑘_𝑓)𝜃^𝑇_𝑖(𝑘_𝑓)}. (95) Due to the additive uncertainties considered in (87), it can be assumed that the average parameter variation is in the range of the same order of magnitude of the fault system. Consequently, it is reasonable to set 𝛿≈1, which denotes the worst case of this assumption. Some numerical examples are also given in [9] to show the sensitivity of 𝛿 versus the mean value of the tracking error to verify the effectiveness of this rule of thumb for resetting the covariance matrix of the parameter estimate.

To improve the parameter estimation for unanticipated faults, 𝑅_𝑖 in (90) needs to be reset by a moving window with the residual𝑅_𝑖≈1𝑘_𝑓−𝑘_𝑖+1^𝑘𝑓_{𝑘=𝑘𝑖}𝜀_𝑖(𝑘)𝜀^𝑇_𝑖(𝑘), (96) where usually 1<(𝑘_𝑓−𝑘_𝑖)<5. Similarly, 𝑄_𝑖 in (89c) should be reset by substituting (95) and (96) into (89b) and then (89a) to obtain𝑄_𝑖≈[𝑀_𝑖(𝑘_𝑓)𝜀_𝑖(𝑘_𝑓)][𝑀_𝑖(𝑘_𝑓)𝜀_𝑖(𝑘_𝑓)]^𝑇. (97)

In the STC scheme, the residual estimate is updated for every sampling time. Consequently, it is convenient to use it as the information of fault detection. Therefore, the time instant 𝑘_𝑓𝑇 of the fault occurrence could be detected by utilizing the ratio of the moving window with the residual 𝑅_𝑖 in (96) and the average norm of the innovation vectors as 𝑅_𝑖𝑅_𝑓𝑖>𝛾_𝑖, (98) where 𝑅_𝑓𝑖=(1/𝑘_𝑓)∑^𝑘𝑓_𝑘=1𝜀_𝑖(𝑘)𝜀^𝑇_𝑖(𝑘) and 𝛾_𝑖 is a preset threshold.

The following summarizes the FTC using the modified STC methodology with the fault detection and covariance matrices resetting.(1)Apply the conventional STC scheme to the healthy system until it tracks the prespecified trajectory well.(2)Switch the conventional RELS estimation algorithm (12) to the modified Kalman filter estimation algorithm in (89a)–(89e) with the initializations of 𝑄_𝑖, 𝑅_𝑖, and 𝑆_𝑖(0) given by (90)–(92).(3)Apply the modified STC scheme to the controlled stochastic system.(4)Execute the fault detection by means of (98) with a preset threshold.(5)Whenever a fault is detected, reset 𝑄_𝑖, 𝑅_𝑖, and 𝑆_𝑖(𝑘_𝑓) given by (95)–(97).(6)Go to step (3) and repeat for the modified STC process.

6. Illustrative Examples

6.1. Parameters Estimation by RELS Method for a Two-Input-Two-Output System

Assume that the two-input-two-output system is unknown, and choose an appropriate two-input-two-output NARMAX model for self-tuning control. Then, identify the parameter of this model through RELS method based on the initial parameter 𝜃(0) of the model obtained by OKID method. Notice that this NARMAX model could approximate other two-input-two-output systems, not just only for the specific example.

Consider the nonlinear two-input-two-output NARMAX model as follows:𝑦₁(𝑘)=𝑦₁(𝑘−1)−0.5𝑦₁(𝑘−2)+𝑦₁(𝑘−2)𝑦₂(𝑘−2)+40𝑢₁(𝑘−1)−2𝑢₂(𝑘−1)+6𝑢²₁(𝑘−2)+𝑒₁𝑦(𝑘−1),₂(𝑘)=−0.4𝑦₂(𝑘−1)−0.5𝑦₁(𝑘−2)−0.2𝑦₁(𝑘−2)𝑦₂(𝑘−2)cos(𝑦₂(𝑘−1))−20𝑢₂(𝑘−1)−8𝑢²₂(𝑘−2)+𝑒₂(𝑘−1), (99) where 𝑦₁(𝑘) and 𝑦₂(𝑘) are outputs, 𝑢₁(𝑘) and 𝑢₂(𝑘) are inputs, and 𝑒₁(𝑘) and 𝑒₂(𝑘) are noises which are zero-mean Gaussian sequences with variances 𝜎²_𝑒1=𝜎²_𝑒2=0.01. Because the system is unknown, the input data and the measurement data of outputs are demanded to estimate the parameter of the selected NARMAX model. The sampling period 𝑇 is 0.1 second. The identification results of ARMAX and NARMAX models through RELS method are shown in Figures 5 and 6, respectively.

(a)

(b)

Figure 5 

The trajectories of the output 1 for ARMAX and NARMAX models.

(a)

(b)

Figure 6 

The trajectories of the output 2 for ARMAX and NARMAX models.

Compare the parameter identification of RELS method based on the intuitive initial parameter 𝜃(0)=[𝐼₂0_2×9𝐼₂0_2×100_2×12]^𝑇 with the OKID-obtained initial parameter 𝜃(0) as follows.

(1) The Intuitive Initial Parameter 𝜃(0) = [𝐼₂0_2×9𝐼₂0_2×100_2×12]^𝑇
Set directly the initial parameter of NARMAX model to be [𝐼₂0_2×9𝐼₂0_2×100_2×12]^𝑇, and then use RELS method given by (12) to identify this model. The results of identification are shown in Figures 7 and 8. Figure 8 shows that the convergence of all different parameters during the process of system identification does not reach a steady state throughout the whole process. Nevertheless, it can be significantly improved by the proposed approach to be shown in Figure 10.

(a)

(b)

(c)

(d)

Figure 7 

The trajectories of outputs using RELS method without OKID.

Figure 8 

Convergence of parameter estimation of NARMAX model using RELS method without OKID.

(2) The Initial Parameter 𝜽(𝟎) Obtained by OKID Method
Assume that the system is excited by the zero mean white noise input 𝑢(𝑡)=[𝑢₁(𝑡)𝑢₂(𝑡)] with the covariance cov(𝑢)=diag{0.1,0.1}, and the sampling period 𝑇 is 0.1 second. Then the system realization estimates {𝐺_𝑜,𝐻_𝑜,𝐶_𝑜} and observer gain 𝐾_𝑜 are obtained by OKID method. Based on the optimal linearization method, one can obtain the initial parameter 𝜃(0) which is close to the convergent value of 𝜃(𝑘). After the initial parameter 𝜃(0) is obtained, RELS method is applied to identify this NARMAX model. The results of identification are shown in Figures 9 and 10. Figures 5 and 6 show that the system identification ability of the proposed NARMAX model-based RELS is better than the traditional ARMAX model-based RELS. From Figure 7 to Figure 10, they show the shorter estimation time and the better system performances through RELS method with OKID.

(a)

(b)

(c)

(d)

Figure 9 

The trajectories of outputs using RELS method with OKID.

Figure 10 

Convergence of parameter estimation of NARMAX model using RELS method with OKID.

6.2. Active Fault Tolerance Using NARMAX Model-Based State-Space Self-Tuning Control with PWM

The system is a nonlinear relative motion dynamics between leader and follower vehicles for multiple spacecraft formation flying (MSFF). In the following model, usual Newton dynamics assumption is made that each spacecraft is evolved as a point-mass in free space. We further assume that the leader is in a circular orbit around the Earth with the constant angular velocity 𝜔.

A schematic drawing of the considered MSFF system is depicted in Figure 11, where we make the following considerations: (i) the inertial coordinate system {𝑋,𝑌,𝑍} is attached to the center of the Earth; (ii) 𝑅_𝐿∈ℜ³ denotes the position from the origin of {𝑋,𝑌,𝑍} to the leader; (iii) the coordinate system {𝑥_𝐿,𝑦_𝐿,𝑧_𝐿} is attached to the leader with the 𝑥_𝐿 axis pointing in the opposite direction as the tangential velocity, the 𝑦_𝐿 axis pointing along the direction of 𝑅_𝐿, and 𝑧_𝐿 axis being mutually perpendicular to the 𝑥_𝐿𝑦_𝐿-plane; overall, {𝑥_𝐿,𝑦_𝐿,𝑧_𝐿} forms a right-hand coordinate frame; (iv) 𝜌=[𝑥,𝑦,𝑧]^𝑇∈ℜ³ denotes the relative position from the origin of {𝑥_𝐿,𝑦_𝐿,𝑧_𝐿} to the follower.

(a)

(b)

Figure 11 

Schematic representations of the MSFF system.

The nonlinear position dynamics of the leader and follower spacecraft with respect to {𝑋,𝑌,𝑍} are given by𝑚_𝐿𝑅̈_𝐿+𝑚_𝐿(𝑀+𝑚_𝐿𝑅)𝐺^𝐿‖𝑅_𝐿‖³+𝐹_𝐿=𝑢_𝐿,𝑚_𝐹(𝑅̈_𝐿+𝜌̈)+𝑚_𝐹(𝑀+𝑚_𝐹𝑅)𝐺^𝐿+𝜌‖𝑅_𝐿+𝜌‖³+𝐹_𝐹=𝑢_𝐹, (100) where 𝑚_𝐿,𝑚_𝐹∈ℜ_>0 are, respectively, the masses of the leader and follower; 𝐹_𝐿,𝐹_𝐹∈ℜ³ are the disturbance forces, respectively; 𝑢_𝐿=[𝑢_𝐿𝑥𝑢_𝐿𝑦𝑢_𝐿𝑧]^𝑇 and 𝑢_𝐹=[𝑢_𝐹𝑥𝑢_𝐹𝑦𝑢_𝐹𝑧]^𝑇 are the actual control inputs of the leader and follower, respectively; 𝑀=5.974×10²⁴ kg is the mass of the Earth; 𝐺=6.673×10⁻¹¹ N-m²/kg² is the universal gravity constant.

Note that the angular velocity is {0,0,𝜔}, the relative acceleration is given by 𝜌̈={𝑥̈−2𝜔𝑦̇−𝜔²𝑥,𝑦̈+2𝜔𝑥̇−𝜔²𝑦,𝑧}̈, and 𝑅_𝐿={0,‖𝑅_𝐿‖,0} is constant in coordinates {𝑥_𝐿,𝑦_𝐿,𝑧_𝐿}. Reflecting M≫𝑚_𝐿,𝑚_𝐹 and calculating some algebraic manipulations on two differential equations are obtain the dynamic equation describing the position of the follower relative to the leader in {𝑥_𝐿,𝑦_𝐿,𝑧_𝐿}, and 𝐹=𝐹_𝐹−(𝑚_𝐹/𝑚_𝐿)𝐹_𝐿=[𝐹_𝑥𝐹_𝑦𝐹_𝑧]^𝑇 is the total constant disturbance. Choosing the states variable 𝑥=[𝜌^𝑇𝜌̇^𝑇]^𝑇=[𝑥₁𝑥₂𝑥₃𝑥₄𝑥₅𝑥₆]^𝑇 and the control 𝑢=𝑢_𝑓, the vector field of the concerned system is represented as follows:̇=𝑥𝑥(𝑡)=𝑓(𝑥(𝑡))+𝐵𝑢(𝑡)+𝑑(𝑡)₄𝑥(𝑡)₅𝑥(𝑡)₆𝜔(𝑡)₂𝑥₁(𝑡)+2𝜔𝑥₅(𝑡)−𝑀𝐺𝑥₁(𝑡)‖𝑅_𝐿+𝜌‖³𝜔²𝑥₂(𝑡)−2𝜔𝑥₄(𝑡)−𝑀𝐺(𝑥₂(𝑡)+‖𝑅_𝐿‖)‖𝑅_𝐿+𝜌‖³−𝑀𝐺𝑥₃(𝑡)‖𝑅_𝐿+𝜌‖³+0₃1𝑚_𝐹𝐼₃000−𝑢(𝑡)+1𝑚_𝐿𝑢_𝐿𝑥−1𝑚_𝐹𝐹_𝑥𝑀𝐺‖𝑅_𝐿‖²−1𝑚_𝐿𝑢_𝐿𝑦−1𝑚_𝐹𝐹_𝑦−1𝑚_𝐿𝑢_𝐿𝑧−1𝑚_𝐹𝐹_𝑧,𝑦(𝑡)=𝑦₁(𝑡)𝑦₂(𝑡)𝑦₃(𝑡)^𝑇=𝐶𝑥(𝑡),𝐶=𝐼₃0₃. (101a)(101b)

The ability for a follower to fly in formation may be limited, in part, by the placement error of the on-off thruster nozzles. A fractional placement error makes the thrust in each axis inefficient and even creates a situation where a spacecraft comes unacceptably close to another. A nonideally placed 𝑥_𝐿-axis thrust denoted by 𝑢_{𝐹𝑥Nonideal} is illustrated in Figure 12, which results in harmful influence on the other axes thrust actions. Putting the relationship 𝑢=𝑢_𝐹Nonideal into (101a) alters input matrix 𝐵=𝐵Θ whereΘ=sin(𝜙_𝑥)cos(𝜃_𝑥)sin(𝜙_𝑦)cos(𝜃_𝑦)sin(𝜙_𝑧)cos(𝜃_𝑧)sin(𝜙_𝑥)sin(𝜃_𝑥)sin(𝜙_𝑦)sin(𝜃_𝑦)sin(𝜙_𝑧)sin(𝜃_𝑧)cos(𝜙_𝑥)cos(𝜙_𝑦)cos(𝜙_𝑧). (102)

Figure 12 

Nonideal placement of 𝑢𝑓𝑥 thruster nozzle.

(a)

(b)

(c)

Figure 13 

Output responses of Scenario 1 without fault-tolerant control.

(a)

(b)

(c)

Figure 14 

Output responses of Scenario 1 with fault-tolerant control.

(a)

(b)

(c)

Figure 15 

PWM control inputs of Scenario 1 with fault-tolerant control.

Consider the MSFF system with parameters as follows: 𝑚_𝐿 = 1550 kg, 𝑚_𝐹 = 410 kg, 𝑅_𝐿 = [04.224×10⁷0]^𝑇 m, 𝜔 = 7.272×10⁻⁵ rad/s, 𝐹 = 10⁻⁵ ×[−1.0256.248−2.415]^𝑇 N, and 𝑢_𝐿 = 0_3×1 N; the relative position and velocity are initialized to 𝜌(0) = [00−200]^𝑇 m and 𝜌̇(0) = 0_3×1 m/s, respectively; the desired relative trajectory is selected as𝑟(𝑡)=100cos(4𝜔𝑡)(1−𝑒^−0.05𝑡3)100sin(4𝜔𝑡)(1−𝑒^−0.05𝑡3)0. (103) Besides, we further introduce the nonideal thruster nozzle placements on all axes by assuming 𝜙_𝑥=𝜙_𝑦=𝜃_𝑦=0.5𝜋, 𝜃_𝑥=0,𝜙_𝑧=0.05𝜋, and 𝜃_𝑧=0.1𝜋 rad. Let the sampling period 𝑇 be 0.05 hour. Because the nonlinear system is unknown for designers, we have only information about inputs and outputs. First, choose an appropriate three-input-three-output NARMAX model for the self-tuning control. The STC equips with an active fault-tolerant controller using the modified estimation algorithm to cope with the unanticipated system faults and input failures. Notice that the initial parameter estimates, 𝜃_𝑖(0)𝑖=1,2,3, of the chosen NARMAX model are obtained by OKID method. If a fault is detected at 𝑡=𝑘_𝑓, matrices 𝑄_𝑖, 𝑅_𝑖, and 𝑆_𝑖(𝑘_𝑓) with the resetting parameter 𝛿=1 in (95) are automatically reset. The fault detection threshold is 𝛾₁=𝛾₂=3.5 and 𝛾₃=5 in (98). Notice that this NARMAX model could approximate other three-input-three-output systems, not just only for this example.

The proposed three-input-three-output NARMAX model is given by𝑦_𝑖(𝑘)=𝐹_𝑖(𝑌)+𝐵_𝑖(𝑌)𝑈+𝐷_𝑖(𝑌)𝐸+𝜀_𝑖(𝑘),𝑖=1,2,3, (104) where=𝑦𝑌₁𝒜𝑦₂𝒜𝑦₃𝒜𝑦₁ℬ𝑦₂ℬ𝑦₃ℬ^𝑇𝐹₁,(𝑌)=𝑎₁₁𝑦₁(𝑘−1)+𝑎₁₂𝑦₂(𝑘−1)+𝑎₁₃𝑦₃(𝑘−1)+𝑎₁₄𝑦₁(𝑘−2)+𝑎₁₅𝑦₂(𝑘−2)+𝑎₁₆𝑦₃(𝑘−2)+𝑎₁₇𝑦²₁(𝑘−1)+𝑎₁₈𝑦²₁(𝑘−2)+𝑎₁₉𝑦₁(𝑘−1)𝑦₂(𝑘−1)+𝑎_1,10𝑦₁(𝑘−1)𝑦₃(𝑘−1)+𝑎_1,11𝑦₁(𝑘−2)𝑦₂(𝑘−1)+𝑎_1,12𝑦₁(𝑘−2)𝑦₃(𝑘−1)+𝑎_1,13𝑦₁(𝑘−1)𝑦₂(𝑘−2)+𝑎_1,14𝑦₁(𝑘−1)𝑦₃(𝑘−2)+𝑎_1,15𝑦₁(𝑘−2)𝑦₂(𝑘−2)+𝑎_1,16𝑦₁(𝑘−2)𝑦₃(𝑘−2)+𝑎_1,17𝑦₁(𝑘−1)𝑦₁𝐹(𝑘−2),₂(𝑌)=𝑎₂₁𝑦₁(𝑘−1)+𝑎₂₂𝑦₂(𝑘−1)+𝑎₂₃𝑦₃(𝑘−1)+𝑎₂₄𝑦₁(𝑘−2)+𝑎₂₅𝑦₂(𝑘−2)+𝑎₂₆𝑦₃(𝑘−2)+𝑎₂₇𝑦²₂(𝑘−1)+𝑎₂₈𝑦²₂(𝑘−2)+𝑎₂₉𝑦₂(𝑘−1)𝑦₁(𝑘−1)+𝑎_2,10𝑦₂(𝑘−1)𝑦₃(𝑘−1)+𝑎_2,11𝑦₂(𝑘−2)𝑦₁(𝑘−1)+𝑎_2,12𝑦₂(𝑘−2)𝑦₃(𝑘−1)+𝑎_2,13𝑦₂(𝑘−1)𝑦₁(𝑘−2)+𝑎_2,14𝑦₂(𝑘−1)𝑦₃(𝑘−2)+𝑎_2,15𝑦₂(𝑘−2)𝑦₁(𝑘−2)+𝑎_2,16𝑦₂(𝑘−2)𝑦₃(𝑘−2)+𝑎_2,17𝑦₂(𝑘−1)𝑦₂𝐹(𝑘−2),₃(𝑌)=𝑎₃₁𝑦₁(𝑘−1)+𝑎₃₂𝑦₂(𝑘−1)+𝑎₃₃𝑦₃(𝑘−1)+𝑎₃₄𝑦₁(𝑘−2)+𝑎₃₅𝑦₂(𝑘−2)+𝑎₃₆𝑦₃(𝑘−2)+𝑎₃₇𝑦²₃(𝑘−1)+𝑎₃₈𝑦²₃(𝑘−2)+𝑎₃₉𝑦₃(𝑘−1)𝑦₁(𝑘−1)+𝑎_3,10𝑦₃(𝑘−1)𝑦₂(𝑘−1)+𝑎_3,11𝑦₃(𝑘−2)𝑦₁(𝑘−1)+𝑎_3,12𝑦₃(𝑘−2)𝑦₂(𝑘−1)+𝑎_3,13𝑦₃(𝑘−1)𝑦₁(𝑘−2)+𝑎_3,14𝑦₃(𝑘−1)𝑦₂(𝑘−2)+𝑎_3,15𝑦₃(𝑘−2)𝑦₁(𝑘−2)+𝑎_3,16𝑦₃(𝑘−2)𝑦₂(𝑘−2)+𝑎_3,17𝑦₃(𝑘−1)𝑦₃𝐵(𝑘−2),_𝑇1(=𝑏𝑌)₁₁+𝑏₁₇𝑦₁(𝑘−1)+𝑏_1,10𝑦₁𝑏(𝑘−2)₁₂+𝑏₁₈𝑦₁(𝑘−1)+𝑏_1,11𝑦₁𝑏(𝑘−2)₁₃+𝑏₁₉𝑦₁(𝑘−1)+𝑏_1,12𝑦₁𝑏(𝑘−2)₁₄+𝑏_1,13𝑦₁(𝑘−1)+𝑏_1,16𝑦₁𝑏(𝑘−2)₁₅+𝑏_1,14𝑦₁(𝑘−1)+𝑏_1,17𝑦₁𝑏(𝑘−2)₁₆+𝑏_1,15𝑦₁(𝑘−1)+𝑏_1,18𝑦₁𝐵(𝑘−2)_𝑇2,(=𝑏𝑌)₂₁+𝑏₂₇𝑦₂(𝑘−1)+𝑏_2,10𝑦₂𝑏(𝑘−2)₂₂+𝑏₂₈𝑦₂(𝑘−1)+𝑏_2,11𝑦₂𝑏(𝑘−2)₂₃+𝑏₂₉𝑦₂(𝑘−1)+𝑏_2,12𝑦₂𝑏(𝑘−2)₂₄+𝑏_2,13𝑦₂(𝑘−1)+𝑏_2,16𝑦₂𝑏(𝑘−2)₂₅+𝑏_2,14𝑦₂(𝑘−1)+𝑏_2,17𝑦₂𝑏(𝑘−2)₂₆+𝑏_2,15𝑦₂(𝑘−1)+𝑏_2,18𝑦₂𝐵(𝑘−2)_𝑇3,(=𝑏𝑌)₃₁+𝑏₃₇𝑦₃(𝑘−1)+𝑏_3,10𝑦₃𝑏(𝑘−2)₃₂+𝑏₃₈𝑦₃(𝑘−1)+𝑏_3,11𝑦₃𝑏(𝑘−2)₃₃+𝑏₃₉𝑦₃(𝑘−1)+𝑏_3,12𝑦₃𝑏(𝑘−2)₃₄+𝑏_3,13𝑦₃(𝑘−1)+𝑏_3,16𝑦₃𝑏(𝑘−2)₃₅+𝑏_3,14𝑦₃(𝑘−1)+𝑏_3,17𝑦₃𝑏(𝑘−2)₃₆+𝑏_3,15𝑦₃(𝑘−1)+𝑏_3,18𝑦₃𝐷(𝑘−2)_𝑇1,(=𝑑𝑌)₁₁+𝑑₁₇𝑦₁(𝑘−1)+𝑑_1,10𝑦₁𝑑(𝑘−2)₁₂+𝑑₁₈𝑦₁(𝑘−1)+𝑑_1,11𝑦₁𝑑(𝑘−2)₁₃+𝑑₁₉𝑦₁(𝑘−1)+𝑑_1,12𝑦₁𝑑(𝑘−2)₁₄+𝑑_1,13𝑦₁(𝑘−1)+𝑑_1,16𝑦₁𝑑(𝑘−2)₁₅+𝑑_1,14𝑦₁(𝑘−1)+𝑑_1,17𝑦₁𝑑(𝑘−2)₁₆+𝑑_1,15𝑦₁(𝑘−1)+𝑑_1,18𝑦₁𝐷(𝑘−2)_𝑇2,(=𝑑𝑌)₂₁+𝑑₂₇𝑦₂(𝑘−1)+𝑑_2,10𝑦₂𝑑(𝑘−2)₂₂+𝑑₂₈𝑦₂(𝑘−1)+𝑑_2,11𝑦₂𝑑(𝑘−2)₂₃+𝑑₂₉𝑦₂(𝑘−1)+𝑑_2,12𝑦₂𝑑(𝑘−2)₂₄+𝑑_2,13𝑦₂(𝑘−1)+𝑑_2,16𝑦₂𝑑(𝑘−2)₂₅+𝑑_2,14𝑦₂(𝑘−1)+𝑑_2,17𝑦₂𝑑(𝑘−2)₂₆+𝑑_2,15𝑦₂(𝑘−1)+𝑑_2,18𝑦₂𝐷(𝑘−2)_𝑇3,(=𝑑𝑌)₃₁+𝑑₃₇𝑦₃(𝑘−1)+𝑑_3,10𝑦₃𝑑(𝑘−2)₃₂+𝑑₃₈𝑦₃(𝑘−1)+𝑑_3,11𝑦₃𝑑(𝑘−2)₃₃+𝑑₃₉𝑦₃(𝑘−1)+𝑑_3,12𝑦₃𝑑(𝑘−2)₃₄+𝑑_3,13𝑦₃(𝑘−1)+𝑑_3,16𝑦₃𝑑(𝑘−2)₃₅+𝑑_3,14𝑦₃(𝑘−1)+𝑑_3,17𝑦₃𝑑(𝑘−2)₃₆+𝑑_3,15𝑦₃(𝑘−1)+𝑑_3,18𝑦₃=𝑢(𝑘−2),𝑈₁𝒜𝑢₂𝒜𝑢₃𝒜𝑢₁ℬ𝑢₂ℬ𝑢₃ℬ^𝑇=𝜀,𝐸₁𝒜𝜀₂𝒜𝜀₃𝒜𝜀₁ℬ𝜀₂ℬ𝜀₃ℬ^𝑇. (105)(27)(28)(29)(30)(31)(32)(33)(34)(35)(36)(37)(38)(39)(40)(41)(42)(43)(44)(45)(46)(47)(48)(49)(50)(51)(52)(53)(54)(55)(56)(57)(58)(59)(60)(61)(62)(63)(64)(65)(66)(67)(68)(69)(70)

Simplify the model to form a linear-in-the-parameter expression as follows:𝑦_𝑖(𝑘)=𝜃^𝑇_𝑖𝜙_𝑖(𝑘)+𝜀_𝑖(𝑘),𝑖=1,2,3, (106) where𝜃_𝑖=𝑎_𝑖1𝑎_𝑖2⋯𝑎_𝑖,17𝑏_𝑖1𝑏_𝑖2⋯𝑏_𝑖,18𝑑_𝑖1𝑑_𝑖2⋯𝑑_𝑖,18^𝑇∈ℜ⁵³.. (107)𝜙_𝑖(𝑘) is composed of the corresponding variables associated with each element of the parameter vector 𝜃_𝑖. 𝜉_𝑖=53 is the number of 𝜃_𝑖. Estimate the parameter 𝜃_𝑖 by using the standard recursive extended-least-squares algorithm, where 𝜆₀=𝜆(0)=0.9,𝑆_𝑖(0)=10×𝐼_𝜉𝑖, and the initial parameter 𝜃_𝑖(0) could be obtained via OKID method.

By the optimal linearization method, the optimal linear models of 𝐹_𝑖(𝑌) at sampling time 𝑘𝑇 are gotten as𝐹_𝑖(𝑌)=𝐴_𝑖𝑌,𝑖=1,2,3, (108) where𝐴_𝑖=𝑎_𝑖1𝑎_𝑖2𝑎_𝑖3𝑎_𝑖4𝑎_𝑖5𝑎_𝑖6𝑎₁₁,𝑎₁₂𝑎₁₃𝑎₁₄𝑎₁₅𝑎₁₆=𝑎₁₁𝑎₁₂𝑎₁₃𝑎₁₄𝑎₁₅𝑎₁₆+2𝑎₁₇𝑎₁₉𝑎_1,10𝑎_1,17𝑎_1,13𝑎_1,14𝑎₁₉00𝑎_1,11𝑎00_1,1000𝑎_1,12𝑎00_1,17𝑎_1,11𝑎_1,122𝑎₁₈𝑎_1,15𝑎_1,16𝑎_1,1300𝑎_1,15𝑎00_1,1400𝑎_1,16𝑌00^𝑜𝑝−𝜂₁𝑌_𝑜𝑝𝜂₁,=1‖𝑌‖²₂×{𝑦₁(𝑘−1)×[𝑎₁₇𝑦₁(𝑘−1)+𝑎₁₉𝑦₂(𝑘−1)+𝑎_1,10𝑦₃(𝑘−1)+𝑎_1,13𝑦₂(𝑘−2)+𝑎_1,14𝑦₃(𝑘−2)+𝑎_1,17𝑦₁(𝑘−2)]+𝑦₁(𝑘−2)×[𝑎₁₈𝑦₁(𝑘−2)+𝑎_1,11𝑦₂(𝑘−1)+𝑎_1,12𝑦₃(𝑘−1)+𝑎_1,15𝑦₂(𝑘−2)+𝑎_1,16𝑦₃(𝑘−2)]}|_𝑌=𝑌_𝑜𝑝𝑎₂₁,𝑎₂₂𝑎₂₃𝑎₂₄𝑎₂₅𝑎₂₆=𝑎₂₁𝑎₂₂𝑎₂₃𝑎₂₄𝑎₂₅𝑎₂₆+0𝑎₂₉00𝑎_2,110𝑎₂₉2𝑎₂₇𝑎_2,10𝑎_2,13𝑎_2,17𝑎_2,140𝑎_2,1000𝑎_2,1200𝑎_2,1300𝑎_2,150𝑎_2,11𝑎_2,17𝑎_2,12𝑎_2,152𝑎₂₈𝑎_2,160𝑎_2,1400𝑎_2,160𝑌_𝑜𝑝−𝜂₂𝑌_𝑜𝑝𝜂₂,=1‖𝑌‖²₂×{𝑦₂(𝑘−1)×[𝑎₂₉𝑦₁(𝑘−1)+𝑎₂₇𝑦₂(𝑘−1)+𝑎_2,10𝑦₃(𝑘−1)+𝑎_2,17𝑦₂(𝑘−2)+𝑎_2,14𝑦₃(𝑘−2)+𝑎_2,13𝑦₁(𝑘−2)]+𝑦₂(𝑘−2)×[𝑎_2,15𝑦₁(𝑘−2)+𝑎_2,11𝑦₁(𝑘−1)+𝑎_2,12𝑦₃(𝑘−1)+𝑎₂₈𝑦₂(𝑘−2)+𝑎_2,16𝑦₃(𝑘−2)]}|_𝑌=𝑌_𝑜𝑝𝑎₃₁,𝑎₃₂𝑎₃₃𝑎₃₄𝑎₃₅𝑎₃₆=𝑎₃₁𝑎₃₂𝑎₃₃𝑎₃₄𝑎₃₅𝑎₃₆+00𝑎₃₉00𝑎_3,1100𝑎_3,1000𝑎_3,12𝑎₃₉𝑎_3,102𝑎₃₇𝑎_3,13𝑎_3,14𝑎_3,1700𝑎_3,1300𝑎_3,1500𝑎_3,1400𝑎_3,16𝑎_3,11𝑎_3,12𝑎_3,17𝑎_3,15𝑎_3,162𝑎₃₈𝑌_𝑜𝑝−𝜂₃𝑌_𝑜𝑝𝜂₃,=1‖𝑌‖²₂×{𝑦₃(𝑘−1)×[𝑎₃₉𝑦₁(𝑘−1)+𝑎_3,10𝑦₂(𝑘−1)+𝑎₃₇𝑦₃(𝑘−1)+𝑎_3,14𝑦₂(𝑘−2)+𝑎_3,17𝑦₃(𝑘−2)+𝑎_3,13𝑦₁(𝑘−2)]+𝑦₃(𝑘−2)×[𝑎_3,15𝑦₁(𝑘−2)+𝑎_3,11𝑦₁(𝑘−1)+𝑎_3,12𝑦₂(𝑘−1)+𝑎_3,16𝑦₂(𝑘−2)+𝑎₃₈𝑦₃(𝑘−2)]}|_𝑌=𝑌_𝑜𝑝. (109)(71)(72)(73)(74)(75)(76)(77)(78)(79)(80)(81)(82)(83)(84)(85)(86)(87)(88)(89)(90)(91)(92)(93)(94)(95)(96)(97)(98)(99)(100)(101)(102)(103)𝐺_𝑜(𝑘) and 𝐻_𝑜(𝑘) for the associated state-space innovation form are determined to be𝐺_𝑜(𝑘)=−𝐺_𝑜1(𝑘)𝐼₃−𝐺_𝑜2(𝑘)0₃,𝐻_𝑜𝐻(𝑘)=^𝑜1𝐻(𝑘)_𝑜2𝐶(𝑘)_𝑜=𝐼₃0₃, (110) where𝐺_𝑜1−𝑎(𝑘)=¹¹−𝑎₁₂−𝑎₁₃−𝑎₂₁−𝑎₂₂−𝑎₂₃−𝑎₃₁−𝑎₃₂−𝑎₃₃,𝐺_𝑜2−𝑎(𝑘)=¹⁴−𝑎₁₅−𝑎₁₆−𝑎₂₄−𝑎₂₅−𝑎₂₆−𝑎₃₄−𝑎₃₅−𝑎₃₆,𝐻_𝑜1𝑏(𝑘)=¹¹+𝑏₁₇𝑦₁𝒜+𝑏_1,10𝑦₁(𝑘−2)𝑏₁₂+𝑏₁₈𝑦₁𝒜+𝑏_1,11𝑦₁(𝑘−2)𝑏₁₃+𝑏₁₉𝑦₁𝒜+𝑏_1,12𝑦₁𝑏(𝑘−2)₂₁+𝑏₂₇𝑦₂𝒜+𝑏_2,10𝑦₂(𝑘−2)𝑏₂₂+𝑏₂₈𝑦₂𝒜+𝑏_2,11𝑦₂(𝑘−2)𝑏₂₃+𝑏₂₉𝑦₂𝒜+𝑏_2,12𝑦₂𝑏(𝑘−2₃₁+𝑏₃₇𝑦₃𝒜+𝑏_3,10𝑦₃(𝑘−2)𝑏₃₂+𝑏₃₈𝑦₃𝒜+𝑏_3,11𝑦₃(𝑘−2)𝑏₃₃+𝑏₃₉𝑦₃𝒜+𝑏_3,12𝑦₃|(𝑘−2)_𝑌=𝑌_𝑜𝑝,𝐻_𝑜2𝑏(𝑘)=¹⁴+𝑏_1,13𝑦₁𝒜+𝑏_1,16𝑦₁(𝑘−2)𝑏₁₅+𝑏_1,14𝑦₁𝒜+𝑏_1,17𝑦₁(𝑘−2)𝑏₁₆+𝑏_1,15𝑦₁𝒜+𝑏_1,18𝑦₁𝑏₂₄+𝑏_2,13𝑦₂𝒜+𝑏_2,16𝑦₂(𝑘−2)𝑏₂₅+𝑏_2,14𝑦₂𝒜+𝑏_2,17𝑦₂(𝑘−2)𝑏₂₆+𝑏_2,15𝑦₂𝒜+𝑏_2,18𝑦₂𝑏₃₄+𝑏_3,13𝑦₃𝒜+𝑏_3,16𝑦₃(𝑘−2)𝑏₃₅+𝑏_3,14𝑦₃𝒜+𝑏_3,17𝑦₃(𝑘−2)𝑏₃₆+𝑏_3,15𝑦₃𝒜+𝑏_3,18𝑦₃|_𝑌=𝑌_𝑜𝑝. (𝑘−2)(𝑘−2)(𝑘−2)(111)The weighting matrices are 𝑄=10⁸×diag{7000,7000,6}, 𝑅=𝐼₃ in (5) and 𝑄_𝑜=10⁵×𝐼₆,𝑅_𝑜=𝐼₃ in (72a)-(72b). The FTC with an abrupt input fault and a gradual input fault is considered.

Scenario (an abrupt input fault). In the beginning, one could obtain the initial parameter estimates, 𝜃_𝑖(0) for 𝑖=1,2,3, of the chosen NARMAX model by OKID method, and use the conventional RELS to identify this model, and then apply the FTC with the modified STC methodology at the 5th hour to fit Assumption 3. At the 11th hour, the Input 1 and Input 2 are assumed to be abruptly added to 2 and 4 times of its function, which fits Assumption 1, respectively. Figure 13 shows that the tracking fails when the system does not equip with the FTC. On the contrary, the tracking of the controlled system with the FTC is relatively better, as shown in Figure 14. The PWM control inputs gotten from the PAM control inputs are shown in Figure 15, which fits Assumption 2.

Scenario (a gradual input fault). In the case, FTC with modified STC methodology is applied at the 5th hour. After the 11th hour, the Input 1 increases by the function 2.5+(𝑡/26)(1−2𝑒^{−6(𝑡−11)}). From the corresponding simulation results, one can see that Figure 16 shows the worse results when the system does not involve the FTC method. On the contrary, the tracking of the system with the FTC is relatively better, as shown in Figure 17. The PWM control inputs gotten from the PAM control inputs are shown in Figure 18.

(a)

(b)

(c)

Figure 16 

Output responses of Scenario 2 without fault-tolerant control.

(a)

(b)

(c)

Figure 17 

Output responses of Scenario 2 with fault-tolerant control.

(a)

(b)

(c)

Figure 18 

PWM control inputs of Scenario 2 with fault-tolerant control.

7. Conclusion

An active fault-tolerant pulse-width-modulated tracker using the NARMAX model-based state-space self-tuner is proposed for continuous-time nonlinear stochastic systems with unknown system parameters and noises. The initial parameter 𝜃_𝑖(0) of the parameter identification of RELS algorithm is determined by using OKID method, rather than based on the intuitive direct-assigned. Hence, a better performance observer can be obtained. Moreover, combining the fault-tolerant control with STC, the modified state-space self-tuning control methodology based on the NARMAX model and OKID method can quickly make an appropriate reaction to the variations of system parameters when the abrupt input fault and/or the gradual input fault occur.

Acknowledgments

This work was supported by the National Science Council of Republic of China under contract NSC96-2221-E-006-292-MY3, NSC98-2221-E-006-159-MY3, the U.S. Army Research Office under Grant W911NF-06-1-0507, the National Science Foundation under Grant NSF 0717860, and the research contract 1440234.