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

Authors: Chu-Tong Wang; Jason S. H. Tsai; Chia-Wei Chen; You Lin; Shu-Mei Guo; Leang-San Shieh
Journal: Journal of Control Science and Engineering
 (2010)

---

## 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.


---

## Body

## 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 Section2, 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:(1a)x((k+1)T)=Gx(kT)+Hu(kT),(1b)y(kT)=Cx(kT),where x(kT)∈Rn, u(kT)∈Rm, and y(kT)∈Rp are the state variable, control input, and control output, respectively. T stands for the sampling period. For simplifying equations and causing no ambiguity, ∙(kT) will be replaced by ∙(k) in many equations of this paper.It is desired to minimize the following quadratic cost function:(2)Jd=12{∑k=1∞[Cx(k)-Γ(k)]TQd[Cx(k)-Γ(k)]+uT(k)Rdu(k)∑k=1∞[Cx(k)-Γ(k)]TQd[Cx(k)-Γ(k)]},
where Qd≥0, Rd>0 and Γ(k)∈Rp is a reference input. The superscript T 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)u(k)=-Kdx(k)+EdΓ*(k),
where(4)Kd=(HTPH+Rd)-1HTPG,Ed=(HTPH+Rd)-1HT[I-(G-HKd)T]-1CTQd.Γ*(k) = Γ(k+1) for the tracking purpose, and P is the positive definite and symmetric solution to the following Riccati equation:(5)GTPG-P-(GTPH)(HTPH+Rd)-1(HTPG)+CTQdC=0.
If the discrete-time system is time-varying, the digital tracker gains (Kd,Ed) 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(6)ẋ(t)=Ax(t)+∑j=1mB(j)u(j)(kT),
where B(j) is the jth column of B and u(j)(kT) is the jth component of the PAM input vector. The corresponding discrete-time model of the system (6) with a zero-order hold is(7a)x(k+1)=Gx(k)+∑j=1mH(j)u(j)(k),
where(7b)G=eAT,H(j)={(G-I)A-1B(j),ifA-1exists,∑i=0∞Ti+1(i+1)!AiB(j),ifA-1doesnotexist.Here I denotes the identity matrix with the same size as G.Let a PWM controlled sampled-data system be represented as(8)ẋw(t)=Axw(t)+∑j=1mB(j)uw(j)(λ),
where the PWM controller is(9)uw(j)(λ)={0,kT≤λ<kT+τk(j),u̅M(j),kT+τk(j)≤λ<kT+τk(j)+δk(j),0,kT+τk(j)+δk(j)≤λ<kT+T,forj=1,2,…,m.
Here u̅M(j), τk(j), and δk(j) are the amplitude of the jth component of the input vector, the firing delay, and the firing duration in the PWM mode at the kth time step, respectively. The graphical illustration of PAM and PWM inputs is shown in Figure 1. The values of τk(j) and δk(j) can be determined as follows:(10)δk(j)=Tu(j)(kT)u̅M(j),τk(j)=(T-δk(j))2.
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 Figure2.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(u(k),u(k-1),u(k-2),…) and the system output data (y(k),y(k-1),y(k-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 θ̂(k-1) is obtained based on the past information up to the time step k-1. This could be used to yield an estimate of the output at the current time step k. 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 θ̂(k-1) to the new value θ̂(k). Depending on the parameter estimate θ̂(k), appropriate controllers can be designed. The reference input Γ(k) 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 them-input p-output system is given by(11)yi(k)=Fi(y(k-1),…,y(k-ny),u(k-1),…,u(k-nu),ε(k-1),…,ε(k-nε)(k-ny))+εi(k)=∑j=1ξiθijϕij(k)+εi(k)=θiTϕi(k)+εi(k),fori=1,2,…,p.
Here yi(k)∈ℜ is the ith output at time step k (k=1,2,3,…); u(k)∈ℜm is the control input; notation εi(k)∈ℜ denotes the ith component of the prediction error ε(k)∈ℜp. 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 ε(k) will be instead of the modeling error. In addition, ny, nu, and nε are the maximum lag indices of y(k), u(k), and ε(k), respectively. θi is the unknown parameter vector for the ith output. ξi denotes the number of θi. The vector ϕi(k) is the linear or nonlinear function of (y(k-1),…,y(k-ny),u(k-1),…,u(k-nu),ε(k-1),…,ε(k-nε)) related to the ith output. Although Fi(·) may be linear or nonlinear polynomial in terms of the past values of (y(k),u(k),ε(k)), the whole model is nonlinear. ϕi(k) may include any linear or nonlinear variables such as terms raised to an integer of power (e.g., y13(k-5), u32(k-4)), products of present and past inputs (e.g., u(k)u(k-2)), products of past outputs (e.g., y2(k-2)y3(k-4)), or cross-terms (e.g., y4(k-1)u2(k-3)).To determine the unknown parameter vectorθi, the standard RELS algorithm can be applied to the linear-in-the-parameter model (11). This RELS algorithm is given by(12)λ(k)=λ0λ(k-1)+(1-λ0),εi(k)=yi(k)-θ̂iT(k-1)ϕi(k),θi(k)=θi(k-1)+Si(k-1)ϕi(k)λ(k)+ϕiT(k)Si(k-1)ϕi(k)εi(k),Si(k)=1λ(k){Si(k-1)-Si(k-1)ϕi(k)ϕiT(k)Si(k-1)λ(k)+ϕiT(k)Si(k-1)ϕi(k)}.
Here λ(k) 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<λ0<1, while λ(k)=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 λ(k) is a tradeoff between high robustness against noises (large λ(k)) and fast tracking capability (small λ(k)). The prediction error εi(k) is the difference between the measured output yi(k) and the one-step-ahead prediction ŷi(k∣k-1)=θ̂iT(k-1)ϕi(k) of yi(k) made at time step k-1 based on the model corresponding to the estimate θ̂i(k-1).If the prediction error is small, the estimateθ̂i(k-1) is “good” and should not be modified very much. The term Si(k) is the parameter estimation error covariance matrix with Si(0)=γiIξi, where Iξi is called the identity matrix of order ξi, and γi is a sufficiently large positive number in order that the parameter estimate can quickly jump away from θ̂i(0).Owing to various kinds of combinations of (y(k),u(k),ε(k)) for Fi(·), 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.


## 2.1. Digital Linear Suboptimal Tracker Design

Consider the underlying linear discrete-time system:(1a)x((k+1)T)=Gx(kT)+Hu(kT),(1b)y(kT)=Cx(kT),where x(kT)∈Rn, u(kT)∈Rm, and y(kT)∈Rp are the state variable, control input, and control output, respectively. T stands for the sampling period. For simplifying equations and causing no ambiguity, ∙(kT) will be replaced by ∙(k) in many equations of this paper.It is desired to minimize the following quadratic cost function:(2)Jd=12{∑k=1∞[Cx(k)-Γ(k)]TQd[Cx(k)-Γ(k)]+uT(k)Rdu(k)∑k=1∞[Cx(k)-Γ(k)]TQd[Cx(k)-Γ(k)]},
where Qd≥0, Rd>0 and Γ(k)∈Rp is a reference input. The superscript T 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)u(k)=-Kdx(k)+EdΓ*(k),
where(4)Kd=(HTPH+Rd)-1HTPG,Ed=(HTPH+Rd)-1HT[I-(G-HKd)T]-1CTQd.Γ*(k) = Γ(k+1) for the tracking purpose, and P is the positive definite and symmetric solution to the following Riccati equation:(5)GTPG-P-(GTPH)(HTPH+Rd)-1(HTPG)+CTQdC=0.
If the discrete-time system is time-varying, the digital tracker gains (Kd,Ed) 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(6)ẋ(t)=Ax(t)+∑j=1mB(j)u(j)(kT),
where B(j) is the jth column of B and u(j)(kT) is the jth component of the PAM input vector. The corresponding discrete-time model of the system (6) with a zero-order hold is(7a)x(k+1)=Gx(k)+∑j=1mH(j)u(j)(k),
where(7b)G=eAT,H(j)={(G-I)A-1B(j),ifA-1exists,∑i=0∞Ti+1(i+1)!AiB(j),ifA-1doesnotexist.Here I denotes the identity matrix with the same size as G.Let a PWM controlled sampled-data system be represented as(8)ẋw(t)=Axw(t)+∑j=1mB(j)uw(j)(λ),
where the PWM controller is(9)uw(j)(λ)={0,kT≤λ<kT+τk(j),u̅M(j),kT+τk(j)≤λ<kT+τk(j)+δk(j),0,kT+τk(j)+δk(j)≤λ<kT+T,forj=1,2,…,m.
Here u̅M(j), τk(j), and δk(j) are the amplitude of the jth component of the input vector, the firing delay, and the firing duration in the PWM mode at the kth time step, respectively. The graphical illustration of PAM and PWM inputs is shown in Figure 1. The values of τk(j) and δk(j) can be determined as follows:(10)δk(j)=Tu(j)(kT)u̅M(j),τk(j)=(T-δk(j))2.
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 Figure2.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(u(k),u(k-1),u(k-2),…) and the system output data (y(k),y(k-1),y(k-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 θ̂(k-1) is obtained based on the past information up to the time step k-1. This could be used to yield an estimate of the output at the current time step k. 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 θ̂(k-1) to the new value θ̂(k). Depending on the parameter estimate θ̂(k), appropriate controllers can be designed. The reference input Γ(k) 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 them-input p-output system is given by(11)yi(k)=Fi(y(k-1),…,y(k-ny),u(k-1),…,u(k-nu),ε(k-1),…,ε(k-nε)(k-ny))+εi(k)=∑j=1ξiθijϕij(k)+εi(k)=θiTϕi(k)+εi(k),fori=1,2,…,p.
Here yi(k)∈ℜ is the ith output at time step k (k=1,2,3,…); u(k)∈ℜm is the control input; notation εi(k)∈ℜ denotes the ith component of the prediction error ε(k)∈ℜp. 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 ε(k) will be instead of the modeling error. In addition, ny, nu, and nε are the maximum lag indices of y(k), u(k), and ε(k), respectively. θi is the unknown parameter vector for the ith output. ξi denotes the number of θi. The vector ϕi(k) is the linear or nonlinear function of (y(k-1),…,y(k-ny),u(k-1),…,u(k-nu),ε(k-1),…,ε(k-nε)) related to the ith output. Although Fi(·) may be linear or nonlinear polynomial in terms of the past values of (y(k),u(k),ε(k)), the whole model is nonlinear. ϕi(k) may include any linear or nonlinear variables such as terms raised to an integer of power (e.g., y13(k-5), u32(k-4)), products of present and past inputs (e.g., u(k)u(k-2)), products of past outputs (e.g., y2(k-2)y3(k-4)), or cross-terms (e.g., y4(k-1)u2(k-3)).To determine the unknown parameter vectorθi, the standard RELS algorithm can be applied to the linear-in-the-parameter model (11). This RELS algorithm is given by(12)λ(k)=λ0λ(k-1)+(1-λ0),εi(k)=yi(k)-θ̂iT(k-1)ϕi(k),θi(k)=θi(k-1)+Si(k-1)ϕi(k)λ(k)+ϕiT(k)Si(k-1)ϕi(k)εi(k),Si(k)=1λ(k){Si(k-1)-Si(k-1)ϕi(k)ϕiT(k)Si(k-1)λ(k)+ϕiT(k)Si(k-1)ϕi(k)}.
Here λ(k) 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<λ0<1, while λ(k)=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 λ(k) is a tradeoff between high robustness against noises (large λ(k)) and fast tracking capability (small λ(k)). The prediction error εi(k) is the difference between the measured output yi(k) and the one-step-ahead prediction ŷi(k∣k-1)=θ̂iT(k-1)ϕi(k) of yi(k) made at time step k-1 based on the model corresponding to the estimate θ̂i(k-1).If the prediction error is small, the estimateθ̂i(k-1) is “good” and should not be modified very much. The term Si(k) is the parameter estimation error covariance matrix with Si(0)=γiIξi, where Iξi is called the identity matrix of order ξi, and γi is a sufficiently large positive number in order that the parameter estimate can quickly jump away from θ̂i(0).Owing to various kinds of combinations of (y(k),u(k),ε(k)) for Fi(·), 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θ̂i(k) 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 θ̂i(0) of the NARMAX model will affect the convergent speed of RELS process. In order to get a suitable initial parameter θ̂i(0) to shorten the transient process of RELS, one could apply the observer/Kalman filter identification (OKID) to evaluate it here.The regression vectorϕi(k) in (11) is composed of(13)[yT𝒜yTℬ⋯yT(k-ny)uT𝒜uTℬ⋯uT(k-nu)εT𝒜εTℬ⋯εT(k-nε)],where 𝒜 denotes (k-1) and ℬ denotes (k-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 u(k) could make the output y(k) of the unknown system track exactly the prespecified reference Γ(k).


### 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(14)x(k+1)=Gx(k)+Hu(k)+w(k),y(k)=Cx(k)+v(k),
where G∈ℜn×n, H∈ℜn×m, and C∈ℜp×n are system, input, and output matrices, respectively; x(k)∈ℜn, u(k)∈ℜm, and y(k)∈ℜp are state, input, and output vectors, respectively; w∈ℜn and v∈ℜp are assumed to be stationary white noise processes with zero mean values and the covariance matrix:(15)E̅([w(k)v(k)][wT(l)vT(l)])=[QSSTR]δ(k-l),fork,l=0,1,2,…,
where E̅(·) is termed the expected value of a random sequence {·}. Q and R are nonnegative and positive definite symmetric matrices for all k, respectively. Due to the assumption that both white noise processes are stationary, the covariance matrix of the process noise, Q, and the covariance matrix of the measurement noise, R, 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(16)Θ=[(CGr-1)T(CGr-2)T⋯(CG)TCT]T
is of full rank.Note that the observability index ofΘ is r=n/p, if it is an integer (otherwise, it is undefined). This constraint means that the Kronecker indices of the system are all such integer r that satisfies n=rp. When this system (14) is block observable, it can be transformed into the block observable companion form as follows:(17)xo(k+1)=Goxo(k)+Hou(k)+wo(k),y(k)=Coxo(k)+vo(k),
where(18)To1=Θ-1[Ip0p⋯0p]T∈ℜn×p,To=[Gr-1To1Gr-2To1⋯GTo1To1]∈ℜn×n,xo(k)=[xo1T(k)xo2T(k)⋯xorT(k)]T=To-1x(k),Go=To-1GTo=[-Go1Ip0p⋯0p-Go20pIp⋯0p⋮⋮⋮⋱⋮-Go(r-1)0p0p⋯Ip-Gor0p0p⋯0p],Ho=To-1H=[Ho1THo2T⋯HorT]T,Co=CTo=[Co1Co2⋯Cor]=[Ip0p⋯0p]∈ℜp×rp,wo(k)=To-1w(k),vo(k)=v(k),E̅{[wo(k)vo(k)][woT(l)voT(l)]}=[To-1Q(To-1)TTo-1S(To-1S)TR]δ(k-l)≜[QoSoSoTRo]δ(k-l),
in which 0p and Ip stand for the p×p null matrix and identity matrix, respectively.The system (17) can be represented by a state-space innovation model [13] as follows:(19)x̂o(k+1∣k)=Gox̂o(k∣k-1)+Hou(k)+Ko(k)eo(k),eo(k)=y(k)-Cox̂o(k∣k-1).
Here the Kalman gain Ko(k) can be computed by the following algorithm [14]:(20)Ko(k)=[GoJo(k)CoT+So][CoJo(k)CoT+Ro]-1,Jo(k+1)=[Go-Ko(k)Co]Jo(k)[Go-Ko(k)Co]T+Ko(k)RoKoT(k)-SoKoT(k)-Ko(k)SoT+Qo=E̅{[xo(k+1)-x̂o(k+1∣k)]T[xo(k+1)-x̂o(k+1∣k)]×[xo(k+1)-x̂o(k+1∣k)]T},Jo(0)=E̅{[xo(0)-x̂o(0)][xo(0)-x̂o(0)]T},
in which x̂o(k∣k-1) is the optimal estimate of xo(k) given by the measurement data y(i) and control input data u(i) for i=0,1,…,k-1; eo(k) is called the innovation process with zero mean and the covariance matrix E̅{eo(k)eoT(k)}=CoJo(k)CoT+Ro.If the pair(Go,Co) is detectable and the pair (Go-SoRo-1Co,Q̅o) is completely stabilizable for any Q̅o satisfying with Q̅oQ̅oT=Qo-SoRo-1SoT [15], then Jo(k)→Jo(∞), where Jo(∞) is the stationary error covariance matrix, so that Ko(k)→Ko(∞) (the stationary Kalman gain) as k→∞. Furthermore, the eigenvalues of Go-Ko(∞)Co are all inside the unit circle centered at the origin. Let z-1 be the backward time-shift operator and set Ko=[Ko1T,Ko2T,…,KorT]T∈ℜn×p. Then, the input-output relationship of the steady-state innovation representation of (19) can be rewritten as (21a)y(k)=Co[In-Goz-1]-1[Hoz-1u(k)+Ko(∞)z-1eo(k)]+eo(k)=[Gl(z-1)]-1[Hl(z-1)]u(k)+[Gl(z-1)]-1[Dl(z-1)]eo(k),
or(21b)Gl(z-1)y(k)=Hl(z-1)u(k)+Dl(z-1)eo(k),where(22)Gl(z-1)=Ip+Go1z-1+Go2z-2+⋯+Gorz-r,Hl(z-1)=Ho1z-1+Ho2z-2+⋯+Horz-r,Dl(z-1)=Ip+Do1z-1+Do2z-2+⋯+Dorz-r,Doi=Goi+Koi,i=1,2,…,r.Note that all the zeros ofdet(Dl(z-1)) 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 Go, Ho, and Co 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 Ko(k). Thus, the state xo(k) can be optimally estimated by using the state-space innovation model (19). When the matrices Go and Ho are unknown, and the covariance matrix is unavailable, the following model in the block observable form(23)x̅o(k+1∣k)=Ĝo(k∣k-1)x̅o(k∣k-1)+Ĥo(k)u(k)+K̂o(k)e̅o(k),e̅o(k)=y(k)-Cox̅o(k∣k-1),
can be used in conjunction with the RELS algorithm to find the Kalman gain estimate K̂o(k) and the state estimate x̅o(k). Here Ĝo(k), Ĥo(k), and K̂o(k) contain parameter estimates Ĝoi(k), Ĥoi(k), and K̂oi(k) for i = 1,2,…,r. When all these parameter estimates converge to the true values, the state estimate x̅o(k∣k-1) of the state xo(k) and the innovation process e̅o(k) converge to the optimal state estimate x̂o(k∣k-1) and innovation process eo(k), respectively. The linear regression model of (23) can be given in the following form:(24a)y(k)=θ̂T(k-1)ϕ(k)+ε(k),
where(24b)θ̂(k)=[Ĝo1(k)Ĝo2(k)⋯Ĝor(k)Ĥo1(k)Ĥo2(k)⋯Ĥor(k)D̂o1(k)D̂o2(k)⋯D̂or(k)]T,ϕ(k)=[-yT𝒜-yTℬ⋯-yT(k-r)uT𝒜uTℬ⋯uT(k-r)εT𝒜εTℬ⋯εT(k-r)]T.It is important to note that as long as the matrix(25)Gc=Go-KoCo=[-Do1Ip0p⋯0p-Do20pIp⋯0p⋮⋮⋮⋱⋮-Do(r-1)0p0p⋯Ip-Dor0p0p⋯0p]
is asymptotically stable, the boundedness of the noise sequences implies that the estimation error will always be bounded. Whenever Doi(k)=0p, for i=1,2,…,r, 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 F, the system (1a) can be rewritten as(26a)x(k+1)=Gx(k)+Hu(k)+Fy(k)-Fy(k)=(G+FC)x(k)+Hu(k)-Fy(k)=G̅x(k)+H̅v(k),
where(26b)G̅=G+FC,(26c)H̅=[H-F],(26d)v(k)=[u(k)y(k)].For a nonzero initial condition,x(0)≠0, the approach should be taken. For this purpose, (26a)–(26d) are extended into the following form:(27)x(k+1)=G̅x(k)+H̅v(k),x(k+2)=G̅x(k+1)+H̅v(k+1)=G̅2x(k)+G̅H̅v(k)+H̅v(k+1),⋮x(k+q)=G̅x(k+q-1)+H̅v(k+q-1)=G̅qx(k)+G̅q-1H̅v(k)+G̅q-2H̅v(k+1)+⋯+G̅H̅v(k+q-2)+H̅v(k+q-1).
Then, the output equation (1b) at time step k+q is given by(28)y(k+q)=Cx(k+q)=CG̅qx(k)+CG̅q-1H̅v(k)+⋯+CG̅H̅v(k+q-2)+CH̅v(k+q-1).
Hence, these output equations, for k=0,1,…,N-q with N>q(m+p+1)-1, can be expressed as(29a)y̅=CG̅qx̅+ϑv̅,
where(29b)y̅=[y(q)y(q+1)⋯y(N)],x̅=[x(0)x(1)⋯x(N-q)],ϑ=[CH̅CG̅H̅⋯CG̅(q-1)H̅],v̅=[v(q-1)v(q)⋯v(N-1)v(q-2)v(q-1)⋯v(N-2)⋮⋮⋱⋮v(0)v(1)⋯v(N-q)].Note that the first term on the right-hand side of (29a) represents the effect of the preceding time steps, where G̅q is sufficiently small and all the states in x are bounded; so (29a) can be approximated by neglecting the first term on the right-hand side, such that(30)y̅≅ϑv̅.
It has the following least-squares solution:(31)ϑ=y̅v̅+,
where v̅+ is the pseudoinverse of the matrix v̅ with v̅v̅+=Iq. It clearly shows that once information of inputs and outputs is obtained, the observer Markov parameters Y̅k=CG̅k-1H̅, for k=1,2,3,…, can be known.The Hankel matrixHa(k-1) [3] is defined by(32)Ha(k-1)≜[Y̅kY̅k+1⋯Y̅k+β-1Y̅k+1Y̅k+2⋯Y̅k+β⋮⋮⋱⋮Y̅k+α-1Y̅k+α⋯Y̅k+α+β-2],k=1,2,3,…,
where α and β are sufficiently large positive integers. Furthermore, Ha(k-1) can be easily decomposed into(33)Ha(k-1)=ΘHaG̅k-1ΦHa,
where ΘHa is the observability matrix defined by(34a)ΘHa=[CT(CG̅)T(CG̅2)T⋯(CG̅α-1)T]T,
and ΦHa is the controllability matrix given by(34b)ΦHa=[H̅G̅H̅G̅2H̅⋯G̅β-1H̅].From (33), one knows that (35a)Ha(0)=ΘHaΦHa,fork=1,(35b)Ha(1)=ΘHaG̅ΦHa,fork=2.When the observer Markov parameters are determined by using (31), the eigensystem realization algorithm (ERA) method is used to obtain the system realization estimates {Ĝ,Ĥ,Ĉ} and the observer gain F through the singular value decomposition (SVD) of the Hankel matrix. The ERA processes the SVD factorization of the block data matrix Ha(k-1), started for k=1, that is,(36)Ha(0)=VΣST.
Here V and S are orthogonal matrices and Σ is a block diagonal matrix of the following form:(37)Σ=[Σñ000],
where Σñ=diag{σ1σ2⋯σnminσnmin+1⋯σñ} contains monotonically nonincreasing entries σ1≥σ2≥⋯≥σnmin≥σnmin+1≥⋯σñ≥0. Here some singular values σnmin+1,…,σñ are relatively small (σnmin+1≪σnmin) 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(38)Σnmin=[σ10⋯000σ2⋯00⋮⋮⋱⋮⋮00⋯σnmin-1000⋯0σnmin].
Hence, (36) can be rewritten as(39)Ha(0)=VnminΣnminSnminT=(VnminΣnmin1/2)(Σnmin1/2SnminT)=ΘHaΦHa,
where Vnmin and Snmin are composed of the first nmin columns of V and S, respectively.Similarly, applying (39) to (35b) gets(40)Ha(1)=ΘHaG̅ΦHa=VnminΣnmin1/2G̅Σnmin1/2SnminT.
Hence, one obtains(41)G̅=Σnmin-1/2VnminTHa(1)SnminΣnmin-1/2.
In other words, the reduced model of order nmin after deleting singular values σnmin+1,…,σñ 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 {Ĝ,Ĥ,Ĉ} and the observer gain Fcan be found out as follows:(42)Ĥ=the  firstmcolumns  ofΣnmin1/2SnminT,-F=the(m+1)th  column  to  the(m+p)th  columnofΣnmin1/2SnminT,Ĉ=the  firstprows  ofVnminΣnmin1/2,Ĝ=Σnmin-1/2VnminTHa(1)SnminΣnmin-1/2-FĈ.
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(43)x(k+1)=f(x(k))+q(x(k))u(k),
where f:ℜn→ℜn and q:ℜn→ℜn×m are nonlinear with continuous partial derivatives with respect to the state variable at all steps k, where x(k)∈ℜn is the state vector at time index k, and u(k)∈ℜm is the control input vector at time index k. It is desired to have an exact local linearization system at an operating state of interest, say xop(k); namely,(44)x(k+1)=G(k)x(k)+H(k)u(k),
where G(k) and H(k) are matrices of appropriate dimensions depending on xop(k). The linearization of the nonlinear system (43) is commonly represented by the truncated Taylor expansion as(45)x(k+1)=f(xeq(k))+q(xeq(k))ueq(k)+A(k)[x(k)-xeq(k)]+B(k)[u(k)-ueq(k)],
where {xeq(k),ueq(k)} is an equilibrium point, A(k)=∇xf(xeq(k)), and B(k)=∇xq(xeq(k)). ∇xf(xeq(k)):ℜn→ℜn×n is the gradient of f(x(k)) (evaluated at x(k)=xeq(k)) with respect to x(k). One can then represent (45) in the following form:(46)x(k+1)=A(k)x(k)+B(k)u(k)-[A(k)xeq(k)+B(k)ueq(k)].
Clearly, (46) is an affine rather than linear model due to the generally nonvanished term A(k)xeq(k)+B(k)ueq(k). One exception is the trivial case where the equilibrium is zero, (xeq(k),ueq(k))=(0,0), which cannot, however, be ensured throughout a nonlinear control process.Suppose thatxop(k)≠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 x(k) and u(k), so that it can well approximate the dynamical behaviors of the system in the vicinity of the operating state xop(k). In other words, one has(47)f(x(k))+q(x(k))u(k)≈G(k)x(k)+H(k)u(k),for anyu(k),f(xop(k))+q(xop(k))u(k)=G(k)xop(k)+H(k)u(k),for anyu(k).
Since the designed control input u(k) is arbitrary, the following equality must be held:(48)q(xop(k))=H(k),
so that (47) becomes quite simple; that is,(49)f(x(k))≈G(k)x(k),f(xop(k))=G(k)xop(k).
To satisfy these, let giT denote the ith row of the matrix G(k). Then (49) can be represented as(50)fi(x(k))≈giTx(k),i=1,2,…,n,(51)fi(xop(k))=giTxop(k),i=1,2,…,n,
where fi:ℜn→ℜ is the ith component of f. Next, expanding the left-hand side of (50) about xop(k) and neglecting the second- and higher-order terms, one has(52)fi(xop(k))+[∇xfi(xop(k))]T(x(k)-xop(k))≈giTx(k).
Now, using (51), one can rewrite (52) as(53)[∇xfi(xop(k))]T(x(k)-xop(k))≈giT(x(k)-xop(k)),
in which x(k) is arbitrary but should be “close” to xop(k) so that the approximation is good. To determine a constant vector giT, such that it is “as close as possible” to [∇xfi(xop(k))]T and also satisfies giTxop(k)=fi(xop(k)), one may consider the following constrained minimization problem:(54)mingiΨ:=12∥∇xfi(xop(k))-gi∥22
subject to(55)giTxop(k)=fi(xop(k)).
Notice that this is a convex constrained optimization problem; therefore, the first-order necessary condition for a minimum of Ψ is also sufficient, which is (56)∇giΨ+η∇gi{giTxop(k)-fi(xop(k))}=0,
where η is the Lagrange multiplier. It follows from (56) that (57)gi-∇xfi(xop(k))+ηxop(k)=0.
Left multiplying (57) by xopT(k), one obtains(58)η=xopT(k)∇xfi(xop(k))-fi(xop(k))∥xop(k)∥22,forxop(k)≠0.
Substituting (58) into (57) gives(59a)gi=∇xfi(xop(k))+fi(xop(k))-xopT(k)∇xfi(xop(k))∥xop(k)∥22xop(k),forxop(k)≠0.
It is easily verified that when xop(k)=0, (57) yields (59b)gi=∇xfi(0),forxop(k)=0.The controllability matrixΦ for the nonlinear system (43) at the operating state xop(k) is derived from the optimal linear model (G(k),H(k)) defined by (44), resulting in(60)Φ=[H̃(k)G̃(k)H̃(k)G̃2(k)H̃(k)⋯G̃n-1(k)H̃(k)],
where G̃(k) and H̃(k) are constructed via the following rule: the ith columns of G(k) and H(k) are set to be zero whenever the ith components of xop(k) and of u(k) 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θ̂i(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 u(k) for simply controlling the unknown true nonlinear system. The elements Go, Ho, and Co of the control gains Kd(k) and Ed(k) 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 withm-inputs, p-outputs and, ρ-time-steps could be selected as follows:(61)yi(k)=Fi(Ỹ)+Bi(Ỹ)U+Di(Ỹ)E+εi(k),fori=1,2,…,p,
where(62)y(k-l)=[y1(k-l)y2(k-l)⋯yp(k-l)]T,forl=1,2,…,ρ,Ỹ=[yT(k-1)yT(k-2)⋯yT(k-ρ)]T∈ℜρp×1,Bi(Ỹ)=[bi1(Ỹ),…,bim(Ỹ),bi(m+1)(Ỹ),…,bi(2m)(Ỹ),…,bi,(ρ-1)×m+1(Ỹ),…,bi(ρ×m)(Ỹ)]∈ℜ1×ρm,Di(Ỹ)=[di1(Ỹ),…,dip(Ỹ),di(p+1)(Ỹ),…,di(2p)(Ỹ),…,di,(ρ-1)×p+1(Ỹ),…,di(ρ×p)(Ỹ)]∈ℜ1×ρp,u(k-l)=[u1(k-l)u2(k-l)⋯um(k-l)]T∈ℜm×1,U=[uT(k-1)uT(k-2)⋯uT(k-ρ)]T∈ℜρm×1,ε(k-l)=[ε1(k-l)ε2(k-l)⋯εp(k-l)]T∈ℜp×1,E=[εT(k-1)εT(k-2)⋯εT(k-ρ)]T∈ℜρp×1.
The prediction errors (εi(k-1),εi(k-2),…,εi(0)) all have been identified from the NARMAX model based on the given control inputs (uj(k-1),uj(k-2),…,uj(0)) and the output measurement data (yi(k),yi(k-1),…,yi(0)) for i=1,2,…,p and j=1,2,…,m. This chosen NARMAX model (61) has the advantage that Bi(Ỹ)U and Di(Ỹ)E are linear in all items of u and ε, respectively.To separate linear/nonlinear variables, this NARMAX model (61) can be rewritten as a linear-in-the-parameter expression as follows:(63)yi(k)=Fi(Ỹ)+Bi(Ỹ)U+Di(Ỹ)E+εi(k)=∑j=1ξiθijϕij(k)+εi(k)=θiTϕi(k)+εi(k),
where(64)θi=[θi1θi2⋯θiξi]T,ϕi(k)=[ϕi1(k)ϕi2(k)⋯ϕiξi(k)]T.Estimate the parameter vectorθi by the RELS algorithm to make the output of the model approximate the true system output. Although Fi(Ỹ) only includes output delays, it is a nonlinear function of y. Therefore, Fi(Ỹ) needs to be linearized by the optimal linearization method in order to find the suboptimal controller.Performing the optimal linearization approach onFi(Ỹ) gives(65a)Fi(Ỹ)=ÃiỸ,
where (65b)ÃiT=(∇ỸFi(Ỹ)+Fi(Ỹ)-ỸT∇ỸFi(Ỹ)∥Ỹ∥22Ỹ)|Ỹ=Ỹop≜[ãi1⋯ãipãi(p+1)⋯ãi(2×p)⋯ãi,(ρ-1)×p+1⋯ãi(ρ×p)]T∈ℜρp×1.Substituting (65a) into (63) with(66)B̃i=Bi(Ỹ)|Ỹ=Ỹop≜[b̃i1,…,b̃im,b̃i(m+1),…,b̃i(2m),…,b̃i,(ρ-1)×m+1,…,b̃i(ρ×m)],D̃i=Di(Ỹ)|Ỹ=Ỹop≜[d̃i1,…,d̃ip,d̃i(p+1),…,d̃i(2p),…,d̃i,(ρ-1)×p+1,…,d̃i(ρ×p)]gives
(67)yi(k)-ÃiỸ=B̃iU+D̃iE+εi(k).
Let z-1 denote the backward shift operator. Equation (67) can be rewritten as(68)yi(k)-ãi1z-1y1(k)-⋯-ãipz-1yp(k)-ãi(p+1)z-2y1(k)+⋯-ãi(2p)z-2yp(k)-⋯-ãi,(ρ-1)×p+1z-ρy1(k)-⋯-ãi(ρ×p)z-ρyp(k)=b̃i1z-1u1(k)+⋯+b̃imz-1um(k)+b̃i(m+1)z-2u1(k)+⋯+b̃i(2×m)z-2um(k)+⋯+b̃i,(ρ-1)×m+1z-ρu1(k)+⋯+b̃i(ρ×m)z-ρum(k)+εi(k)+d̃i1z-1ε1(k)+⋯+d̃ipz-1εp(k)+d̃i(p+1)z-2ε1(k)+⋯+d̃i(2×p)z-2εp(k)+⋯+d̃i,(ρ-1)×p+1z-ρε1(k)+⋯+d̃i(ρ×p)z-ρεp(k),fori=1,2,…,p.After combining (68) for i=1,2,…,p, a disturbed output of the plant in the ARMAX model takes the following form:(69a)GL(z-1)y(k)=HL(z-1)u(k)+DL(z-1)ε(k),
where
(69b)GL(z-1)=Ip+Go1z-1+⋯+Goρz-ρ,Gol(l=1,2,…,ρ)∈ℜp×p,HL(z-1)=Ho1z-1+Ho2z-2+⋯+Hoρz-ρ,Hol(l=1,2,…,ρ)∈ℜp×m,DL(z-1)=Ip+Do1z-1+⋯+Doρz-ρ,Dol(l=1,2,…,ρ)∈ℜp×p,Gol=[-ã1,(l-1)×p+1-ã1,(l-1)×p+2⋯-ã1(l×p)-ã2,(l-1)×p+1-ã2,(l-1)×p+2⋯-ã2(l×p)⋮⋮⋱⋮-ãp,(l-1)×p+1-ãp,(l-1)×p+2⋯-ãp(l×p)],Hol=[b̃1,(l-1)×m+1b̃1,(l-1)×m+2⋯b̃1(l×m)b̃2,(l-1)×m+1b̃2,(l-1)×m+2⋯b̃2(l×m)⋮⋮⋱⋮b̃p,(l-1)×m+1b̃p,(l-1)×m+2⋯b̃p(l×m)],Dol=[d̃1,(l-1)×p+1d̃1,(l-1)×p+2⋯d̃1(l×p)d̃2,(l-1)×p+1d̃2,(l-1)×p+2⋯d̃2(l×p)⋮⋮⋱⋮d̃p,(l-1)×p+1d̃p,(l-1)×p+2⋯d̃p(l×p)].An alternative representation of the ARMAX model (69a) and (69b) is given by(70)y(k)=GL-1(z-1)HL(z-1)u(k)+GL-1(z-1)DL(z-1)ε(k),
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 GL-1(z-1), which represents the effects of control inputs and noises. Once GL-1(z-1) has been specified to characterize the dynamics of the plant, the error vector model GL-1(z-1)DL(z-1)ε(k) presents an adjustable moving average process of DL(z-1)ε(k) 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]:(71a)x̂o(k+1)=Gox̂o(k)+Hou(k)+Ko(k)eo(k),(71b)eo(k)=y(k)-Cox̂o(k),
where(71c)Go=[-Go1Ip0p⋯0p-Go20pIp⋯0p⋮⋮⋮⋱⋮-Go(ρ-1)0p0p⋯Ip-Goρ0p0p⋯0p],(71d)Ho=[Ho1THo2T⋯HoρT]T,(71e)Co=[Ip0p⋯0p],(71f)Ko=[Ko1Ko2⋯Koρ],Kol=Dol-Golforl=1,2,…,ρ,(71g)x̂o(k)=[x̂o1T(k)x̂o2T(k)⋯x̂oρT(k)]T,x̂ol(k)∈ℜpforl=1,2,…,ρ.x̂o(k) is the estimate of the system state xo(k) in the observer coordinates. In this process, the matrices Gol, Hol, and Dol are calculated as Gol(k), Hol(k), and Dol(k) at each sampling time kT, respectively.However, the zeros ofDL(z-1) 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 Ko(k) may not all lie in the unit circle centered at the origin, either. Instead of the Kalman gain Ko(k), one could design the digital estimator gain Lo(k) to replace Ko(k). The digital estimator gain Lo(k) 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 (GoT(k),CoT) is stabilizable, there exists a digital estimator gain Lo(k) so that the closed-loop estimator matrix Go(k)-Lo(k)Co is asymptotically stable. The estimator gain is(72a)LoT(k)=(CoPo(k)CoT+Ro)-1CoPo(k)GoT(k).
Here Po(k) is the solution to the following Riccati equation:(72b)Go(k)Po(k)GoT(k)-Po(k)-(Go(k)Po(k)CoT)×(CoPo(k)CoT+Ro)-1(CoPo(k)GoT(k))+Qo=0,in which Qo≥0 and Ro>0 with appropriate dimensions. Then the state-space innovation form (71a) becomes(73)x̂o(k+1)=Go(k)x̂o(k)+Ho(k)u(k)+Lo(k)eo(k).
Thus, one can design the digital controller(74)u(k)=-Kd(k)x̂o(k)+Ed(k)Γ*(k)
by the observer (73).


#### 3.4.2. The Initial Parameters of NARMAX Model Based on OKID

The initial parameter vectorθ̂i(0) of NARMAX model will affect the convergent speed of RELS process. In order to increase the convergent speed of RELS algorithm, θ̂i(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 {Ĝ,Ĥ,Ĉ} and observer gain F. Then, transfer these matrices into the corresponding observer form {Ĝo,Ĥo,Ĉo} and observer gain Ko (71a)–(71g) by(75)Θ=[(ĈĜρ-1)T(ĈĜρ-2)T⋯(ĈĜ)TĈT]T,To1=Θ-1[Ip0p⋯0p]T∈ℜn×p,To=[Ĝρ-1To1Ĝρ-2To1⋯ĜTo1To1]∈ℜn×n,Ĝo=To-1ĜTo∈ℜn×n,Ĥo=To-1Ĥ=[Ĥo1TĤo2T⋯ĤoρT]T∈ℜn×m,Ĉo=ĈTo=[Ip0p⋯0p]∈ℜp×n,Ko=To-1F=[Ko1Ko2⋯Koρ].
Based on (69a)-(69b) and (71a)–(71g) one could get an ARMAX model as follows:(76)y(k)+Ĝo1y(k-1)+Ĝo2y(k-2)+⋯+Ĝoρy(k-ρ)=Ĥo1u(k-1)+Ĥo2u(k-2)+⋯+Ĥoρu(k-ρ)+D̂o1ε(k-1)+D̂o2ε(k-2)+⋯+D̂oρε(k-ρ),
where D̂ol=Ĝol+Kol,l=1,2,…,ρ.Comparing (76) with (69a) and (69b) obtains the coefficient matrices {Ãi,B̃i,D̃i} of the linearized ARMAX model defined by (65a)-(65b), (66), and (67). Rearranging the right-hand sides of (65b) and (66), respectively, gives(77)Ãi=[∇ỸFi(Ỹop)+Fi(Ỹop)-ỸopT∇ỸFi(Ỹop)∥Ỹop∥22Ỹop]T=[ai1ai2⋯ainFi]F̅i(Ỹop),B̃i=Bi(Ỹ)|Ỹ=Ỹop=[bi1bi2⋯binBi]B̅i(Ỹop),D̃i=Di(Ỹ)|Ỹ=Ỹop=[di1di2⋯dinDi]D̅i(Ỹop),
where i=1,2,…,p. The unknown parameters ail, bim, and din are the undetermined coefficients of Fi(Ỹ),Bi(Ỹ), and Di(Ỹ) for l=1,2,…,nFi, m=1,2,…,nBi, and n=1,2,…,nDi, respectively. F̅i(Ỹop), B̅i(Ỹop), and D̅i(Ỹop) are composed of the corresponding terms with respect to ail, bim, and din, respectively.Note thatξi=nFi+nBi+nDi. 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:(78)[âi1âi2⋯âinFi]=ÃiF̅i+(Ỹop),[b̂i1b̂i2⋯b̂inBi]=B̃iB̅i+(Ỹop),[d̂i1d̂i2⋯d̂inDi]=D̃iD̅i+(Ỹop).
Consequently, the parameter vector estimate(79)θ̂i=[âi1âi2⋯âinFib̂i1b̂i2⋯b̂inBid̂i1d̂i2⋯d̂inDi]T=[ÃiF̅i+(Ỹop)B̃iB̅i+(Ỹop)D̃iD̅i+(Ỹop)]T
can be served as the initial parameter θ̂i(0) of RELS algorithm (12) for STC.


## 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(14)x(k+1)=Gx(k)+Hu(k)+w(k),y(k)=Cx(k)+v(k),
where G∈ℜn×n, H∈ℜn×m, and C∈ℜp×n are system, input, and output matrices, respectively; x(k)∈ℜn, u(k)∈ℜm, and y(k)∈ℜp are state, input, and output vectors, respectively; w∈ℜn and v∈ℜp are assumed to be stationary white noise processes with zero mean values and the covariance matrix:(15)E̅([w(k)v(k)][wT(l)vT(l)])=[QSSTR]δ(k-l),fork,l=0,1,2,…,
where E̅(·) is termed the expected value of a random sequence {·}. Q and R are nonnegative and positive definite symmetric matrices for all k, respectively. Due to the assumption that both white noise processes are stationary, the covariance matrix of the process noise, Q, and the covariance matrix of the measurement noise, R, 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(16)Θ=[(CGr-1)T(CGr-2)T⋯(CG)TCT]T
is of full rank.Note that the observability index ofΘ is r=n/p, if it is an integer (otherwise, it is undefined). This constraint means that the Kronecker indices of the system are all such integer r that satisfies n=rp. When this system (14) is block observable, it can be transformed into the block observable companion form as follows:(17)xo(k+1)=Goxo(k)+Hou(k)+wo(k),y(k)=Coxo(k)+vo(k),
where(18)To1=Θ-1[Ip0p⋯0p]T∈ℜn×p,To=[Gr-1To1Gr-2To1⋯GTo1To1]∈ℜn×n,xo(k)=[xo1T(k)xo2T(k)⋯xorT(k)]T=To-1x(k),Go=To-1GTo=[-Go1Ip0p⋯0p-Go20pIp⋯0p⋮⋮⋮⋱⋮-Go(r-1)0p0p⋯Ip-Gor0p0p⋯0p],Ho=To-1H=[Ho1THo2T⋯HorT]T,Co=CTo=[Co1Co2⋯Cor]=[Ip0p⋯0p]∈ℜp×rp,wo(k)=To-1w(k),vo(k)=v(k),E̅{[wo(k)vo(k)][woT(l)voT(l)]}=[To-1Q(To-1)TTo-1S(To-1S)TR]δ(k-l)≜[QoSoSoTRo]δ(k-l),
in which 0p and Ip stand for the p×p null matrix and identity matrix, respectively.The system (17) can be represented by a state-space innovation model [13] as follows:(19)x̂o(k+1∣k)=Gox̂o(k∣k-1)+Hou(k)+Ko(k)eo(k),eo(k)=y(k)-Cox̂o(k∣k-1).
Here the Kalman gain Ko(k) can be computed by the following algorithm [14]:(20)Ko(k)=[GoJo(k)CoT+So][CoJo(k)CoT+Ro]-1,Jo(k+1)=[Go-Ko(k)Co]Jo(k)[Go-Ko(k)Co]T+Ko(k)RoKoT(k)-SoKoT(k)-Ko(k)SoT+Qo=E̅{[xo(k+1)-x̂o(k+1∣k)]T[xo(k+1)-x̂o(k+1∣k)]×[xo(k+1)-x̂o(k+1∣k)]T},Jo(0)=E̅{[xo(0)-x̂o(0)][xo(0)-x̂o(0)]T},
in which x̂o(k∣k-1) is the optimal estimate of xo(k) given by the measurement data y(i) and control input data u(i) for i=0,1,…,k-1; eo(k) is called the innovation process with zero mean and the covariance matrix E̅{eo(k)eoT(k)}=CoJo(k)CoT+Ro.If the pair(Go,Co) is detectable and the pair (Go-SoRo-1Co,Q̅o) is completely stabilizable for any Q̅o satisfying with Q̅oQ̅oT=Qo-SoRo-1SoT [15], then Jo(k)→Jo(∞), where Jo(∞) is the stationary error covariance matrix, so that Ko(k)→Ko(∞) (the stationary Kalman gain) as k→∞. Furthermore, the eigenvalues of Go-Ko(∞)Co are all inside the unit circle centered at the origin. Let z-1 be the backward time-shift operator and set Ko=[Ko1T,Ko2T,…,KorT]T∈ℜn×p. Then, the input-output relationship of the steady-state innovation representation of (19) can be rewritten as (21a)y(k)=Co[In-Goz-1]-1[Hoz-1u(k)+Ko(∞)z-1eo(k)]+eo(k)=[Gl(z-1)]-1[Hl(z-1)]u(k)+[Gl(z-1)]-1[Dl(z-1)]eo(k),
or(21b)Gl(z-1)y(k)=Hl(z-1)u(k)+Dl(z-1)eo(k),where(22)Gl(z-1)=Ip+Go1z-1+Go2z-2+⋯+Gorz-r,Hl(z-1)=Ho1z-1+Ho2z-2+⋯+Horz-r,Dl(z-1)=Ip+Do1z-1+Do2z-2+⋯+Dorz-r,Doi=Goi+Koi,i=1,2,…,r.Note that all the zeros ofdet(Dl(z-1)) 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 Go, Ho, and Co 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 Ko(k). Thus, the state xo(k) can be optimally estimated by using the state-space innovation model (19). When the matrices Go and Ho are unknown, and the covariance matrix is unavailable, the following model in the block observable form(23)x̅o(k+1∣k)=Ĝo(k∣k-1)x̅o(k∣k-1)+Ĥo(k)u(k)+K̂o(k)e̅o(k),e̅o(k)=y(k)-Cox̅o(k∣k-1),
can be used in conjunction with the RELS algorithm to find the Kalman gain estimate K̂o(k) and the state estimate x̅o(k). Here Ĝo(k), Ĥo(k), and K̂o(k) contain parameter estimates Ĝoi(k), Ĥoi(k), and K̂oi(k) for i = 1,2,…,r. When all these parameter estimates converge to the true values, the state estimate x̅o(k∣k-1) of the state xo(k) and the innovation process e̅o(k) converge to the optimal state estimate x̂o(k∣k-1) and innovation process eo(k), respectively. The linear regression model of (23) can be given in the following form:(24a)y(k)=θ̂T(k-1)ϕ(k)+ε(k),
where(24b)θ̂(k)=[Ĝo1(k)Ĝo2(k)⋯Ĝor(k)Ĥo1(k)Ĥo2(k)⋯Ĥor(k)D̂o1(k)D̂o2(k)⋯D̂or(k)]T,ϕ(k)=[-yT𝒜-yTℬ⋯-yT(k-r)uT𝒜uTℬ⋯uT(k-r)εT𝒜εTℬ⋯εT(k-r)]T.It is important to note that as long as the matrix(25)Gc=Go-KoCo=[-Do1Ip0p⋯0p-Do20pIp⋯0p⋮⋮⋮⋱⋮-Do(r-1)0p0p⋯Ip-Dor0p0p⋯0p]
is asymptotically stable, the boundedness of the noise sequences implies that the estimation error will always be bounded. Whenever Doi(k)=0p, for i=1,2,…,r, 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 F, the system (1a) can be rewritten as(26a)x(k+1)=Gx(k)+Hu(k)+Fy(k)-Fy(k)=(G+FC)x(k)+Hu(k)-Fy(k)=G̅x(k)+H̅v(k),
where(26b)G̅=G+FC,(26c)H̅=[H-F],(26d)v(k)=[u(k)y(k)].For a nonzero initial condition,x(0)≠0, the approach should be taken. For this purpose, (26a)–(26d) are extended into the following form:(27)x(k+1)=G̅x(k)+H̅v(k),x(k+2)=G̅x(k+1)+H̅v(k+1)=G̅2x(k)+G̅H̅v(k)+H̅v(k+1),⋮x(k+q)=G̅x(k+q-1)+H̅v(k+q-1)=G̅qx(k)+G̅q-1H̅v(k)+G̅q-2H̅v(k+1)+⋯+G̅H̅v(k+q-2)+H̅v(k+q-1).
Then, the output equation (1b) at time step k+q is given by(28)y(k+q)=Cx(k+q)=CG̅qx(k)+CG̅q-1H̅v(k)+⋯+CG̅H̅v(k+q-2)+CH̅v(k+q-1).
Hence, these output equations, for k=0,1,…,N-q with N>q(m+p+1)-1, can be expressed as(29a)y̅=CG̅qx̅+ϑv̅,
where(29b)y̅=[y(q)y(q+1)⋯y(N)],x̅=[x(0)x(1)⋯x(N-q)],ϑ=[CH̅CG̅H̅⋯CG̅(q-1)H̅],v̅=[v(q-1)v(q)⋯v(N-1)v(q-2)v(q-1)⋯v(N-2)⋮⋮⋱⋮v(0)v(1)⋯v(N-q)].Note that the first term on the right-hand side of (29a) represents the effect of the preceding time steps, where G̅q is sufficiently small and all the states in x are bounded; so (29a) can be approximated by neglecting the first term on the right-hand side, such that(30)y̅≅ϑv̅.
It has the following least-squares solution:(31)ϑ=y̅v̅+,
where v̅+ is the pseudoinverse of the matrix v̅ with v̅v̅+=Iq. It clearly shows that once information of inputs and outputs is obtained, the observer Markov parameters Y̅k=CG̅k-1H̅, for k=1,2,3,…, can be known.The Hankel matrixHa(k-1) [3] is defined by(32)Ha(k-1)≜[Y̅kY̅k+1⋯Y̅k+β-1Y̅k+1Y̅k+2⋯Y̅k+β⋮⋮⋱⋮Y̅k+α-1Y̅k+α⋯Y̅k+α+β-2],k=1,2,3,…,
where α and β are sufficiently large positive integers. Furthermore, Ha(k-1) can be easily decomposed into(33)Ha(k-1)=ΘHaG̅k-1ΦHa,
where ΘHa is the observability matrix defined by(34a)ΘHa=[CT(CG̅)T(CG̅2)T⋯(CG̅α-1)T]T,
and ΦHa is the controllability matrix given by(34b)ΦHa=[H̅G̅H̅G̅2H̅⋯G̅β-1H̅].From (33), one knows that (35a)Ha(0)=ΘHaΦHa,fork=1,(35b)Ha(1)=ΘHaG̅ΦHa,fork=2.When the observer Markov parameters are determined by using (31), the eigensystem realization algorithm (ERA) method is used to obtain the system realization estimates {Ĝ,Ĥ,Ĉ} and the observer gain F through the singular value decomposition (SVD) of the Hankel matrix. The ERA processes the SVD factorization of the block data matrix Ha(k-1), started for k=1, that is,(36)Ha(0)=VΣST.
Here V and S are orthogonal matrices and Σ is a block diagonal matrix of the following form:(37)Σ=[Σñ000],
where Σñ=diag{σ1σ2⋯σnminσnmin+1⋯σñ} contains monotonically nonincreasing entries σ1≥σ2≥⋯≥σnmin≥σnmin+1≥⋯σñ≥0. Here some singular values σnmin+1,…,σñ are relatively small (σnmin+1≪σnmin) 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(38)Σnmin=[σ10⋯000σ2⋯00⋮⋮⋱⋮⋮00⋯σnmin-1000⋯0σnmin].
Hence, (36) can be rewritten as(39)Ha(0)=VnminΣnminSnminT=(VnminΣnmin1/2)(Σnmin1/2SnminT)=ΘHaΦHa,
where Vnmin and Snmin are composed of the first nmin columns of V and S, respectively.Similarly, applying (39) to (35b) gets(40)Ha(1)=ΘHaG̅ΦHa=VnminΣnmin1/2G̅Σnmin1/2SnminT.
Hence, one obtains(41)G̅=Σnmin-1/2VnminTHa(1)SnminΣnmin-1/2.
In other words, the reduced model of order nmin after deleting singular values σnmin+1,…,σñ 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 {Ĝ,Ĥ,Ĉ} and the observer gain Fcan be found out as follows:(42)Ĥ=the  firstmcolumns  ofΣnmin1/2SnminT,-F=the(m+1)th  column  to  the(m+p)th  columnofΣnmin1/2SnminT,Ĉ=the  firstprows  ofVnminΣnmin1/2,Ĝ=Σnmin-1/2VnminTHa(1)SnminΣnmin-1/2-FĈ.
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(43)x(k+1)=f(x(k))+q(x(k))u(k),
where f:ℜn→ℜn and q:ℜn→ℜn×m are nonlinear with continuous partial derivatives with respect to the state variable at all steps k, where x(k)∈ℜn is the state vector at time index k, and u(k)∈ℜm is the control input vector at time index k. It is desired to have an exact local linearization system at an operating state of interest, say xop(k); namely,(44)x(k+1)=G(k)x(k)+H(k)u(k),
where G(k) and H(k) are matrices of appropriate dimensions depending on xop(k). The linearization of the nonlinear system (43) is commonly represented by the truncated Taylor expansion as(45)x(k+1)=f(xeq(k))+q(xeq(k))ueq(k)+A(k)[x(k)-xeq(k)]+B(k)[u(k)-ueq(k)],
where {xeq(k),ueq(k)} is an equilibrium point, A(k)=∇xf(xeq(k)), and B(k)=∇xq(xeq(k)). ∇xf(xeq(k)):ℜn→ℜn×n is the gradient of f(x(k)) (evaluated at x(k)=xeq(k)) with respect to x(k). One can then represent (45) in the following form:(46)x(k+1)=A(k)x(k)+B(k)u(k)-[A(k)xeq(k)+B(k)ueq(k)].
Clearly, (46) is an affine rather than linear model due to the generally nonvanished term A(k)xeq(k)+B(k)ueq(k). One exception is the trivial case where the equilibrium is zero, (xeq(k),ueq(k))=(0,0), which cannot, however, be ensured throughout a nonlinear control process.Suppose thatxop(k)≠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 x(k) and u(k), so that it can well approximate the dynamical behaviors of the system in the vicinity of the operating state xop(k). In other words, one has(47)f(x(k))+q(x(k))u(k)≈G(k)x(k)+H(k)u(k),for anyu(k),f(xop(k))+q(xop(k))u(k)=G(k)xop(k)+H(k)u(k),for anyu(k).
Since the designed control input u(k) is arbitrary, the following equality must be held:(48)q(xop(k))=H(k),
so that (47) becomes quite simple; that is,(49)f(x(k))≈G(k)x(k),f(xop(k))=G(k)xop(k).
To satisfy these, let giT denote the ith row of the matrix G(k). Then (49) can be represented as(50)fi(x(k))≈giTx(k),i=1,2,…,n,(51)fi(xop(k))=giTxop(k),i=1,2,…,n,
where fi:ℜn→ℜ is the ith component of f. Next, expanding the left-hand side of (50) about xop(k) and neglecting the second- and higher-order terms, one has(52)fi(xop(k))+[∇xfi(xop(k))]T(x(k)-xop(k))≈giTx(k).
Now, using (51), one can rewrite (52) as(53)[∇xfi(xop(k))]T(x(k)-xop(k))≈giT(x(k)-xop(k)),
in which x(k) is arbitrary but should be “close” to xop(k) so that the approximation is good. To determine a constant vector giT, such that it is “as close as possible” to [∇xfi(xop(k))]T and also satisfies giTxop(k)=fi(xop(k)), one may consider the following constrained minimization problem:(54)mingiΨ:=12∥∇xfi(xop(k))-gi∥22
subject to(55)giTxop(k)=fi(xop(k)).
Notice that this is a convex constrained optimization problem; therefore, the first-order necessary condition for a minimum of Ψ is also sufficient, which is (56)∇giΨ+η∇gi{giTxop(k)-fi(xop(k))}=0,
where η is the Lagrange multiplier. It follows from (56) that (57)gi-∇xfi(xop(k))+ηxop(k)=0.
Left multiplying (57) by xopT(k), one obtains(58)η=xopT(k)∇xfi(xop(k))-fi(xop(k))∥xop(k)∥22,forxop(k)≠0.
Substituting (58) into (57) gives(59a)gi=∇xfi(xop(k))+fi(xop(k))-xopT(k)∇xfi(xop(k))∥xop(k)∥22xop(k),forxop(k)≠0.
It is easily verified that when xop(k)=0, (57) yields (59b)gi=∇xfi(0),forxop(k)=0.The controllability matrixΦ for the nonlinear system (43) at the operating state xop(k) is derived from the optimal linear model (G(k),H(k)) defined by (44), resulting in(60)Φ=[H̃(k)G̃(k)H̃(k)G̃2(k)H̃(k)⋯G̃n-1(k)H̃(k)],
where G̃(k) and H̃(k) are constructed via the following rule: the ith columns of G(k) and H(k) are set to be zero whenever the ith components of xop(k) and of u(k) 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θ̂i(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 u(k) for simply controlling the unknown true nonlinear system. The elements Go, Ho, and Co of the control gains Kd(k) and Ed(k) 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 withm-inputs, p-outputs and, ρ-time-steps could be selected as follows:(61)yi(k)=Fi(Ỹ)+Bi(Ỹ)U+Di(Ỹ)E+εi(k),fori=1,2,…,p,
where(62)y(k-l)=[y1(k-l)y2(k-l)⋯yp(k-l)]T,forl=1,2,…,ρ,Ỹ=[yT(k-1)yT(k-2)⋯yT(k-ρ)]T∈ℜρp×1,Bi(Ỹ)=[bi1(Ỹ),…,bim(Ỹ),bi(m+1)(Ỹ),…,bi(2m)(Ỹ),…,bi,(ρ-1)×m+1(Ỹ),…,bi(ρ×m)(Ỹ)]∈ℜ1×ρm,Di(Ỹ)=[di1(Ỹ),…,dip(Ỹ),di(p+1)(Ỹ),…,di(2p)(Ỹ),…,di,(ρ-1)×p+1(Ỹ),…,di(ρ×p)(Ỹ)]∈ℜ1×ρp,u(k-l)=[u1(k-l)u2(k-l)⋯um(k-l)]T∈ℜm×1,U=[uT(k-1)uT(k-2)⋯uT(k-ρ)]T∈ℜρm×1,ε(k-l)=[ε1(k-l)ε2(k-l)⋯εp(k-l)]T∈ℜp×1,E=[εT(k-1)εT(k-2)⋯εT(k-ρ)]T∈ℜρp×1.
The prediction errors (εi(k-1),εi(k-2),…,εi(0)) all have been identified from the NARMAX model based on the given control inputs (uj(k-1),uj(k-2),…,uj(0)) and the output measurement data (yi(k),yi(k-1),…,yi(0)) for i=1,2,…,p and j=1,2,…,m. This chosen NARMAX model (61) has the advantage that Bi(Ỹ)U and Di(Ỹ)E are linear in all items of u and ε, respectively.To separate linear/nonlinear variables, this NARMAX model (61) can be rewritten as a linear-in-the-parameter expression as follows:(63)yi(k)=Fi(Ỹ)+Bi(Ỹ)U+Di(Ỹ)E+εi(k)=∑j=1ξiθijϕij(k)+εi(k)=θiTϕi(k)+εi(k),
where(64)θi=[θi1θi2⋯θiξi]T,ϕi(k)=[ϕi1(k)ϕi2(k)⋯ϕiξi(k)]T.Estimate the parameter vectorθi by the RELS algorithm to make the output of the model approximate the true system output. Although Fi(Ỹ) only includes output delays, it is a nonlinear function of y. Therefore, Fi(Ỹ) needs to be linearized by the optimal linearization method in order to find the suboptimal controller.Performing the optimal linearization approach onFi(Ỹ) gives(65a)Fi(Ỹ)=ÃiỸ,
where (65b)ÃiT=(∇ỸFi(Ỹ)+Fi(Ỹ)-ỸT∇ỸFi(Ỹ)∥Ỹ∥22Ỹ)|Ỹ=Ỹop≜[ãi1⋯ãipãi(p+1)⋯ãi(2×p)⋯ãi,(ρ-1)×p+1⋯ãi(ρ×p)]T∈ℜρp×1.Substituting (65a) into (63) with(66)B̃i=Bi(Ỹ)|Ỹ=Ỹop≜[b̃i1,…,b̃im,b̃i(m+1),…,b̃i(2m),…,b̃i,(ρ-1)×m+1,…,b̃i(ρ×m)],D̃i=Di(Ỹ)|Ỹ=Ỹop≜[d̃i1,…,d̃ip,d̃i(p+1),…,d̃i(2p),…,d̃i,(ρ-1)×p+1,…,d̃i(ρ×p)]gives
(67)yi(k)-ÃiỸ=B̃iU+D̃iE+εi(k).
Let z-1 denote the backward shift operator. Equation (67) can be rewritten as(68)yi(k)-ãi1z-1y1(k)-⋯-ãipz-1yp(k)-ãi(p+1)z-2y1(k)+⋯-ãi(2p)z-2yp(k)-⋯-ãi,(ρ-1)×p+1z-ρy1(k)-⋯-ãi(ρ×p)z-ρyp(k)=b̃i1z-1u1(k)+⋯+b̃imz-1um(k)+b̃i(m+1)z-2u1(k)+⋯+b̃i(2×m)z-2um(k)+⋯+b̃i,(ρ-1)×m+1z-ρu1(k)+⋯+b̃i(ρ×m)z-ρum(k)+εi(k)+d̃i1z-1ε1(k)+⋯+d̃ipz-1εp(k)+d̃i(p+1)z-2ε1(k)+⋯+d̃i(2×p)z-2εp(k)+⋯+d̃i,(ρ-1)×p+1z-ρε1(k)+⋯+d̃i(ρ×p)z-ρεp(k),fori=1,2,…,p.After combining (68) for i=1,2,…,p, a disturbed output of the plant in the ARMAX model takes the following form:(69a)GL(z-1)y(k)=HL(z-1)u(k)+DL(z-1)ε(k),
where
(69b)GL(z-1)=Ip+Go1z-1+⋯+Goρz-ρ,Gol(l=1,2,…,ρ)∈ℜp×p,HL(z-1)=Ho1z-1+Ho2z-2+⋯+Hoρz-ρ,Hol(l=1,2,…,ρ)∈ℜp×m,DL(z-1)=Ip+Do1z-1+⋯+Doρz-ρ,Dol(l=1,2,…,ρ)∈ℜp×p,Gol=[-ã1,(l-1)×p+1-ã1,(l-1)×p+2⋯-ã1(l×p)-ã2,(l-1)×p+1-ã2,(l-1)×p+2⋯-ã2(l×p)⋮⋮⋱⋮-ãp,(l-1)×p+1-ãp,(l-1)×p+2⋯-ãp(l×p)],Hol=[b̃1,(l-1)×m+1b̃1,(l-1)×m+2⋯b̃1(l×m)b̃2,(l-1)×m+1b̃2,(l-1)×m+2⋯b̃2(l×m)⋮⋮⋱⋮b̃p,(l-1)×m+1b̃p,(l-1)×m+2⋯b̃p(l×m)],Dol=[d̃1,(l-1)×p+1d̃1,(l-1)×p+2⋯d̃1(l×p)d̃2,(l-1)×p+1d̃2,(l-1)×p+2⋯d̃2(l×p)⋮⋮⋱⋮d̃p,(l-1)×p+1d̃p,(l-1)×p+2⋯d̃p(l×p)].An alternative representation of the ARMAX model (69a) and (69b) is given by(70)y(k)=GL-1(z-1)HL(z-1)u(k)+GL-1(z-1)DL(z-1)ε(k),
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 GL-1(z-1), which represents the effects of control inputs and noises. Once GL-1(z-1) has been specified to characterize the dynamics of the plant, the error vector model GL-1(z-1)DL(z-1)ε(k) presents an adjustable moving average process of DL(z-1)ε(k) 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]:(71a)x̂o(k+1)=Gox̂o(k)+Hou(k)+Ko(k)eo(k),(71b)eo(k)=y(k)-Cox̂o(k),
where(71c)Go=[-Go1Ip0p⋯0p-Go20pIp⋯0p⋮⋮⋮⋱⋮-Go(ρ-1)0p0p⋯Ip-Goρ0p0p⋯0p],(71d)Ho=[Ho1THo2T⋯HoρT]T,(71e)Co=[Ip0p⋯0p],(71f)Ko=[Ko1Ko2⋯Koρ],Kol=Dol-Golforl=1,2,…,ρ,(71g)x̂o(k)=[x̂o1T(k)x̂o2T(k)⋯x̂oρT(k)]T,x̂ol(k)∈ℜpforl=1,2,…,ρ.x̂o(k) is the estimate of the system state xo(k) in the observer coordinates. In this process, the matrices Gol, Hol, and Dol are calculated as Gol(k), Hol(k), and Dol(k) at each sampling time kT, respectively.However, the zeros ofDL(z-1) 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 Ko(k) may not all lie in the unit circle centered at the origin, either. Instead of the Kalman gain Ko(k), one could design the digital estimator gain Lo(k) to replace Ko(k). The digital estimator gain Lo(k) 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 (GoT(k),CoT) is stabilizable, there exists a digital estimator gain Lo(k) so that the closed-loop estimator matrix Go(k)-Lo(k)Co is asymptotically stable. The estimator gain is(72a)LoT(k)=(CoPo(k)CoT+Ro)-1CoPo(k)GoT(k).
Here Po(k) is the solution to the following Riccati equation:(72b)Go(k)Po(k)GoT(k)-Po(k)-(Go(k)Po(k)CoT)×(CoPo(k)CoT+Ro)-1(CoPo(k)GoT(k))+Qo=0,in which Qo≥0 and Ro>0 with appropriate dimensions. Then the state-space innovation form (71a) becomes(73)x̂o(k+1)=Go(k)x̂o(k)+Ho(k)u(k)+Lo(k)eo(k).
Thus, one can design the digital controller(74)u(k)=-Kd(k)x̂o(k)+Ed(k)Γ*(k)
by the observer (73).


### 3.4.2. The Initial Parameters of NARMAX Model Based on OKID

The initial parameter vectorθ̂i(0) of NARMAX model will affect the convergent speed of RELS process. In order to increase the convergent speed of RELS algorithm, θ̂i(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 {Ĝ,Ĥ,Ĉ} and observer gain F. Then, transfer these matrices into the corresponding observer form {Ĝo,Ĥo,Ĉo} and observer gain Ko (71a)–(71g) by(75)Θ=[(ĈĜρ-1)T(ĈĜρ-2)T⋯(ĈĜ)TĈT]T,To1=Θ-1[Ip0p⋯0p]T∈ℜn×p,To=[Ĝρ-1To1Ĝρ-2To1⋯ĜTo1To1]∈ℜn×n,Ĝo=To-1ĜTo∈ℜn×n,Ĥo=To-1Ĥ=[Ĥo1TĤo2T⋯ĤoρT]T∈ℜn×m,Ĉo=ĈTo=[Ip0p⋯0p]∈ℜp×n,Ko=To-1F=[Ko1Ko2⋯Koρ].
Based on (69a)-(69b) and (71a)–(71g) one could get an ARMAX model as follows:(76)y(k)+Ĝo1y(k-1)+Ĝo2y(k-2)+⋯+Ĝoρy(k-ρ)=Ĥo1u(k-1)+Ĥo2u(k-2)+⋯+Ĥoρu(k-ρ)+D̂o1ε(k-1)+D̂o2ε(k-2)+⋯+D̂oρε(k-ρ),
where D̂ol=Ĝol+Kol,l=1,2,…,ρ.Comparing (76) with (69a) and (69b) obtains the coefficient matrices {Ãi,B̃i,D̃i} of the linearized ARMAX model defined by (65a)-(65b), (66), and (67). Rearranging the right-hand sides of (65b) and (66), respectively, gives(77)Ãi=[∇ỸFi(Ỹop)+Fi(Ỹop)-ỸopT∇ỸFi(Ỹop)∥Ỹop∥22Ỹop]T=[ai1ai2⋯ainFi]F̅i(Ỹop),B̃i=Bi(Ỹ)|Ỹ=Ỹop=[bi1bi2⋯binBi]B̅i(Ỹop),D̃i=Di(Ỹ)|Ỹ=Ỹop=[di1di2⋯dinDi]D̅i(Ỹop),
where i=1,2,…,p. The unknown parameters ail, bim, and din are the undetermined coefficients of Fi(Ỹ),Bi(Ỹ), and Di(Ỹ) for l=1,2,…,nFi, m=1,2,…,nBi, and n=1,2,…,nDi, respectively. F̅i(Ỹop), B̅i(Ỹop), and D̅i(Ỹop) are composed of the corresponding terms with respect to ail, bim, and din, respectively.Note thatξi=nFi+nBi+nDi. 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:(78)[âi1âi2⋯âinFi]=ÃiF̅i+(Ỹop),[b̂i1b̂i2⋯b̂inBi]=B̃iB̅i+(Ỹop),[d̂i1d̂i2⋯d̂inDi]=D̃iD̅i+(Ỹop).
Consequently, the parameter vector estimate(79)θ̂i=[âi1âi2⋯âinFib̂i1b̂i2⋯b̂inBid̂i1d̂i2⋯d̂inDi]T=[ÃiF̅i+(Ỹop)B̃iB̅i+(Ỹop)D̃iD̅i+(Ỹop)]T
can be served as the initial parameter θ̂i(0) of RELS algorithm (12) for STC.


## 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 u(k) for simply controlling the unknown true nonlinear system. The elements Go, Ho, and Co of the control gains Kd(k) and Ed(k) 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 withm-inputs, p-outputs and, ρ-time-steps could be selected as follows:(61)yi(k)=Fi(Ỹ)+Bi(Ỹ)U+Di(Ỹ)E+εi(k),fori=1,2,…,p,
where(62)y(k-l)=[y1(k-l)y2(k-l)⋯yp(k-l)]T,forl=1,2,…,ρ,Ỹ=[yT(k-1)yT(k-2)⋯yT(k-ρ)]T∈ℜρp×1,Bi(Ỹ)=[bi1(Ỹ),…,bim(Ỹ),bi(m+1)(Ỹ),…,bi(2m)(Ỹ),…,bi,(ρ-1)×m+1(Ỹ),…,bi(ρ×m)(Ỹ)]∈ℜ1×ρm,Di(Ỹ)=[di1(Ỹ),…,dip(Ỹ),di(p+1)(Ỹ),…,di(2p)(Ỹ),…,di,(ρ-1)×p+1(Ỹ),…,di(ρ×p)(Ỹ)]∈ℜ1×ρp,u(k-l)=[u1(k-l)u2(k-l)⋯um(k-l)]T∈ℜm×1,U=[uT(k-1)uT(k-2)⋯uT(k-ρ)]T∈ℜρm×1,ε(k-l)=[ε1(k-l)ε2(k-l)⋯εp(k-l)]T∈ℜp×1,E=[εT(k-1)εT(k-2)⋯εT(k-ρ)]T∈ℜρp×1.
The prediction errors (εi(k-1),εi(k-2),…,εi(0)) all have been identified from the NARMAX model based on the given control inputs (uj(k-1),uj(k-2),…,uj(0)) and the output measurement data (yi(k),yi(k-1),…,yi(0)) for i=1,2,…,p and j=1,2,…,m. This chosen NARMAX model (61) has the advantage that Bi(Ỹ)U and Di(Ỹ)E are linear in all items of u and ε, respectively.To separate linear/nonlinear variables, this NARMAX model (61) can be rewritten as a linear-in-the-parameter expression as follows:(63)yi(k)=Fi(Ỹ)+Bi(Ỹ)U+Di(Ỹ)E+εi(k)=∑j=1ξiθijϕij(k)+εi(k)=θiTϕi(k)+εi(k),
where(64)θi=[θi1θi2⋯θiξi]T,ϕi(k)=[ϕi1(k)ϕi2(k)⋯ϕiξi(k)]T.Estimate the parameter vectorθi by the RELS algorithm to make the output of the model approximate the true system output. Although Fi(Ỹ) only includes output delays, it is a nonlinear function of y. Therefore, Fi(Ỹ) needs to be linearized by the optimal linearization method in order to find the suboptimal controller.Performing the optimal linearization approach onFi(Ỹ) gives(65a)Fi(Ỹ)=ÃiỸ,
where (65b)ÃiT=(∇ỸFi(Ỹ)+Fi(Ỹ)-ỸT∇ỸFi(Ỹ)∥Ỹ∥22Ỹ)|Ỹ=Ỹop≜[ãi1⋯ãipãi(p+1)⋯ãi(2×p)⋯ãi,(ρ-1)×p+1⋯ãi(ρ×p)]T∈ℜρp×1.Substituting (65a) into (63) with(66)B̃i=Bi(Ỹ)|Ỹ=Ỹop≜[b̃i1,…,b̃im,b̃i(m+1),…,b̃i(2m),…,b̃i,(ρ-1)×m+1,…,b̃i(ρ×m)],D̃i=Di(Ỹ)|Ỹ=Ỹop≜[d̃i1,…,d̃ip,d̃i(p+1),…,d̃i(2p),…,d̃i,(ρ-1)×p+1,…,d̃i(ρ×p)]gives
(67)yi(k)-ÃiỸ=B̃iU+D̃iE+εi(k).
Let z-1 denote the backward shift operator. Equation (67) can be rewritten as(68)yi(k)-ãi1z-1y1(k)-⋯-ãipz-1yp(k)-ãi(p+1)z-2y1(k)+⋯-ãi(2p)z-2yp(k)-⋯-ãi,(ρ-1)×p+1z-ρy1(k)-⋯-ãi(ρ×p)z-ρyp(k)=b̃i1z-1u1(k)+⋯+b̃imz-1um(k)+b̃i(m+1)z-2u1(k)+⋯+b̃i(2×m)z-2um(k)+⋯+b̃i,(ρ-1)×m+1z-ρu1(k)+⋯+b̃i(ρ×m)z-ρum(k)+εi(k)+d̃i1z-1ε1(k)+⋯+d̃ipz-1εp(k)+d̃i(p+1)z-2ε1(k)+⋯+d̃i(2×p)z-2εp(k)+⋯+d̃i,(ρ-1)×p+1z-ρε1(k)+⋯+d̃i(ρ×p)z-ρεp(k),fori=1,2,…,p.After combining (68) for i=1,2,…,p, a disturbed output of the plant in the ARMAX model takes the following form:(69a)GL(z-1)y(k)=HL(z-1)u(k)+DL(z-1)ε(k),
where
(69b)GL(z-1)=Ip+Go1z-1+⋯+Goρz-ρ,Gol(l=1,2,…,ρ)∈ℜp×p,HL(z-1)=Ho1z-1+Ho2z-2+⋯+Hoρz-ρ,Hol(l=1,2,…,ρ)∈ℜp×m,DL(z-1)=Ip+Do1z-1+⋯+Doρz-ρ,Dol(l=1,2,…,ρ)∈ℜp×p,Gol=[-ã1,(l-1)×p+1-ã1,(l-1)×p+2⋯-ã1(l×p)-ã2,(l-1)×p+1-ã2,(l-1)×p+2⋯-ã2(l×p)⋮⋮⋱⋮-ãp,(l-1)×p+1-ãp,(l-1)×p+2⋯-ãp(l×p)],Hol=[b̃1,(l-1)×m+1b̃1,(l-1)×m+2⋯b̃1(l×m)b̃2,(l-1)×m+1b̃2,(l-1)×m+2⋯b̃2(l×m)⋮⋮⋱⋮b̃p,(l-1)×m+1b̃p,(l-1)×m+2⋯b̃p(l×m)],Dol=[d̃1,(l-1)×p+1d̃1,(l-1)×p+2⋯d̃1(l×p)d̃2,(l-1)×p+1d̃2,(l-1)×p+2⋯d̃2(l×p)⋮⋮⋱⋮d̃p,(l-1)×p+1d̃p,(l-1)×p+2⋯d̃p(l×p)].An alternative representation of the ARMAX model (69a) and (69b) is given by(70)y(k)=GL-1(z-1)HL(z-1)u(k)+GL-1(z-1)DL(z-1)ε(k),
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 GL-1(z-1), which represents the effects of control inputs and noises. Once GL-1(z-1) has been specified to characterize the dynamics of the plant, the error vector model GL-1(z-1)DL(z-1)ε(k) presents an adjustable moving average process of DL(z-1)ε(k) 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]:(71a)x̂o(k+1)=Gox̂o(k)+Hou(k)+Ko(k)eo(k),(71b)eo(k)=y(k)-Cox̂o(k),
where(71c)Go=[-Go1Ip0p⋯0p-Go20pIp⋯0p⋮⋮⋮⋱⋮-Go(ρ-1)0p0p⋯Ip-Goρ0p0p⋯0p],(71d)Ho=[Ho1THo2T⋯HoρT]T,(71e)Co=[Ip0p⋯0p],(71f)Ko=[Ko1Ko2⋯Koρ],Kol=Dol-Golforl=1,2,…,ρ,(71g)x̂o(k)=[x̂o1T(k)x̂o2T(k)⋯x̂oρT(k)]T,x̂ol(k)∈ℜpforl=1,2,…,ρ.x̂o(k) is the estimate of the system state xo(k) in the observer coordinates. In this process, the matrices Gol, Hol, and Dol are calculated as Gol(k), Hol(k), and Dol(k) at each sampling time kT, respectively.However, the zeros ofDL(z-1) 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 Ko(k) may not all lie in the unit circle centered at the origin, either. Instead of the Kalman gain Ko(k), one could design the digital estimator gain Lo(k) to replace Ko(k). The digital estimator gain Lo(k) 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 (GoT(k),CoT) is stabilizable, there exists a digital estimator gain Lo(k) so that the closed-loop estimator matrix Go(k)-Lo(k)Co is asymptotically stable. The estimator gain is(72a)LoT(k)=(CoPo(k)CoT+Ro)-1CoPo(k)GoT(k).
Here Po(k) is the solution to the following Riccati equation:(72b)Go(k)Po(k)GoT(k)-Po(k)-(Go(k)Po(k)CoT)×(CoPo(k)CoT+Ro)-1(CoPo(k)GoT(k))+Qo=0,in which Qo≥0 and Ro>0 with appropriate dimensions. Then the state-space innovation form (71a) becomes(73)x̂o(k+1)=Go(k)x̂o(k)+Ho(k)u(k)+Lo(k)eo(k).
Thus, one can design the digital controller(74)u(k)=-Kd(k)x̂o(k)+Ed(k)Γ*(k)
by the observer (73).


## 3.4.2. The Initial Parameters of NARMAX Model Based on OKID

The initial parameter vectorθ̂i(0) of NARMAX model will affect the convergent speed of RELS process. In order to increase the convergent speed of RELS algorithm, θ̂i(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 {Ĝ,Ĥ,Ĉ} and observer gain F. Then, transfer these matrices into the corresponding observer form {Ĝo,Ĥo,Ĉo} and observer gain Ko (71a)–(71g) by(75)Θ=[(ĈĜρ-1)T(ĈĜρ-2)T⋯(ĈĜ)TĈT]T,To1=Θ-1[Ip0p⋯0p]T∈ℜn×p,To=[Ĝρ-1To1Ĝρ-2To1⋯ĜTo1To1]∈ℜn×n,Ĝo=To-1ĜTo∈ℜn×n,Ĥo=To-1Ĥ=[Ĥo1TĤo2T⋯ĤoρT]T∈ℜn×m,Ĉo=ĈTo=[Ip0p⋯0p]∈ℜp×n,Ko=To-1F=[Ko1Ko2⋯Koρ].
Based on (69a)-(69b) and (71a)–(71g) one could get an ARMAX model as follows:(76)y(k)+Ĝo1y(k-1)+Ĝo2y(k-2)+⋯+Ĝoρy(k-ρ)=Ĥo1u(k-1)+Ĥo2u(k-2)+⋯+Ĥoρu(k-ρ)+D̂o1ε(k-1)+D̂o2ε(k-2)+⋯+D̂oρε(k-ρ),
where D̂ol=Ĝol+Kol,l=1,2,…,ρ.Comparing (76) with (69a) and (69b) obtains the coefficient matrices {Ãi,B̃i,D̃i} of the linearized ARMAX model defined by (65a)-(65b), (66), and (67). Rearranging the right-hand sides of (65b) and (66), respectively, gives(77)Ãi=[∇ỸFi(Ỹop)+Fi(Ỹop)-ỸopT∇ỸFi(Ỹop)∥Ỹop∥22Ỹop]T=[ai1ai2⋯ainFi]F̅i(Ỹop),B̃i=Bi(Ỹ)|Ỹ=Ỹop=[bi1bi2⋯binBi]B̅i(Ỹop),D̃i=Di(Ỹ)|Ỹ=Ỹop=[di1di2⋯dinDi]D̅i(Ỹop),
where i=1,2,…,p. The unknown parameters ail, bim, and din are the undetermined coefficients of Fi(Ỹ),Bi(Ỹ), and Di(Ỹ) for l=1,2,…,nFi, m=1,2,…,nBi, and n=1,2,…,nDi, respectively. F̅i(Ỹop), B̅i(Ỹop), and D̅i(Ỹop) are composed of the corresponding terms with respect to ail, bim, and din, respectively.Note thatξi=nFi+nBi+nDi. 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:(78)[âi1âi2⋯âinFi]=ÃiF̅i+(Ỹop),[b̂i1b̂i2⋯b̂inBi]=B̃iB̅i+(Ỹop),[d̂i1d̂i2⋯d̂inDi]=D̃iD̅i+(Ỹop).
Consequently, the parameter vector estimate(79)θ̂i=[âi1âi2⋯âinFib̂i1b̂i2⋯b̂inBid̂i1d̂i2⋯d̂inDi]T=[ÃiF̅i+(Ỹop)B̃iB̅i+(Ỹop)D̃iD̅i+(Ỹop)]T
can be served as the initial parameter θ̂i(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)ẋ(t)=f(x(t))+q(x(t))u(t)+w(t),y(t)=h(x(t))+v(t),
where f:ℜn→ℜn, q:ℜn→ℜn×m, h:ℜn→ℜp, u(k)∈ℜm is the control input, x(k)∈ℜn is the state vector, y(t)∈ℜp is the measurable output vector, and w(t)∈ℜn and v(t)∈ℜp are uncorrelated white noise processes. It is assumed that f(x(t)), q(x(t)), h(x(t)), w(t), and v(t) are unknown, and x(k) 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 Section3.2 to obtain system and observer-gain Markov parameters of the OKID model; next, use the ERA method to obtain the system realization estimates Ĝ, Ĥ, Ĉ, and observer gain F and then transfer them into observer form Ĝo, Ĥo, Ĉo, and Ko.(ii)
based on the state-space innovation form (71a)–(71g) and the optimal linearization of NARMAX model in (65a)–(69b), the initial parameter θ̂i(0) can be obtained by using (79).Step 2.
Identify the parameter of the NARMAX model at each sampling periodT.(i)
Let the number ofθi be ξi. Also, set Si(0)=γiIξi>0, 0<λ0<1, 0.9<λ(0)<1, x̂o(0), k=1, and θ̂i(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 θ̂i(k) at each sampling interval T by RELS algorithm.Step 3.
Linearize the NARMAX model at each sampling periodT:
Based on the updated parameter estimateθ̂i(k), linearizing the NARMAX model by the optimal linearization methodology and using (69a)-(69b) and (71a)–(71g) yields Go(k), Ho(k), and Co.Step 4.
Estimate the predicted state at each sampling periodT:
Select appropriate{Qo,Ro} in (72b) to have the high-gain property digital estimator gain:(81)LoT(k)=(CoPo(k)CoT+Ro)-1CoPo(k)GoT(k),
where Po(k) is the positive definite and symmetric solution to the Riccati equation (72b). The associated state-space innovation form is given by
(82)x̂o(k∣k-1)=Go(k-1)x̂o(k-1∣k-2)+Ho(k-1)u(k-1)+Lo(k-1)eo(k-1),eo(k-1)=y(k-1)-Cox̂o(k-1∣k-2).
From the above three equations, one can predict the statex̂o(k∣k-1).Step 5.
Generate the digital control input at each sampling periodT.(i)
Select appropriate weighting matrices{Qd,Rd} for the Riccati equation (5) such that the desired digital control law possesses the high-gain property.(ii)
Compute the digital control gainsKd(k) and Ed(k) by the digital control formula (4) as follows:
(83)Kd(k)=[HoT(k)P(k)Ho(k)+Rd]-1HoT(k)P(k)Go(k),Ed(k)=[HoT(k)P(k)Ho(k)+Rd]-1×HoT(k)[I-(Go(k)-Ho(k)Kd)T]-1CoTQd,
where P(k) is the positive definite symmetric solution of the following Riccati equation:
(84)GoT(k)P(k)Go(k)-P(k)+CoTQoCo-GoT(k)P(k)Ho(k)×[HoT(k)P(k)Ho(k)+Ro]-1HoT(k)P(k)Go(k)=0.(iii)
Generate the digital PAM control input at each sampling periodT by using (73) and (74) with Kd(k) and Ed(k) obtained in Step 5(ii)
(85)u(t)=u(kT)=-Kd(kT)x̂o(kT∣kT-T)+Ed(kT)Γ(kT+T),forkT≤t<kT+T.
Note that the reference input Γ(t)=Γ(kT+T), for kT≤t<kT+T, is based on the tracking purpose.(iv)
Calculate the PWM control inputuw(t) by using (9) and (10).(v)
Setk=k+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(86)ẋ(t)=f(x)+q(x)u(t)+∑v=1ςβv(t-τv)[fv(x)+qv(u)]+w(t),
where fv(x)+qv(u) 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 βv(t) as shown in Figure 4 [21], where us(t) denotes the unit-step function of time t. fv(x)+qv(u), βv(t), and τv 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) The gradual failure function. (b) The abrupt failure function.
(a)(b)By the optimal linearization method described in Section3.3, the nonlinear faulty system (86) can be accurately linearized as the following linear time-varying model:(87)ẋ(t)=[G(t)+∑v=1ςβv(t-τv)Gv(t)]x(t)+[H(t)+∑v=1ςβv(t-τv)Hv(t)]u(t)+w(t).
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 τv. 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:(88a)θi(k+1)=θi(k)+w̅i(k),(88b)yi(k)=θiT(k)ϕi(k)+v̅i(k),(88c)E̅(w̅i(k)w̅iT(k))=Qi,(88d)E̅(v̅i(k)v̅iT(k))=Ri,where i=1,2,…,p. {w̅i(k)} and {v̅i(k)} are white Gaussian noise sequences with zero means. Applying the Kalman filter to this model (88a) and (88d) gives the following recursive algorithm:(89a)θ̂i(k)=θ̂i(k-1)+Mi(k)εi(k),(89b)Mi(k)=S̅i(k-1)ϕi(k)Ri+ϕiT(k)S̅i(k-1)ϕi(k),(89c)S̅i(k)=S̅i(k-1)-S̅i(k-1)ϕi(k)ϕiT(k)S̅i(k-1)Ri+ϕiT(k)S̅i(k-1)ϕi(k)+Qi,(89d)εi(k)=yi(k)-θ̂iT(k-1)ϕi(k),
with(89e)S̅i(0)=E̅{[θ̂i(0)-θi(0)][θ̂i(0)-θi(0)]T},where S̅i(k) is the covariance matrix of the parameter estimate error. {εi(k)} 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 ofQi, Ri, and S̅i(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 ofQi, Ri, and S̅i(0) to the following moving window-based statistical quantities: (90)Ri≈1N-k1+1∑k=k1Nεi(k)εiT(k),(91)S̅i(0)≈diag{∑k=k1N[(θ̂i(k)-θi(N))(θ̂i(k)-θi(N))T]diagonal  elements  of1N-k1+1×∑k=k1N[(θ̂i(k)-θi(N))(θ̂i(k)-θi(N))T]},
where k1 is the time index after the estimate θ̂i(k) is in steady state, and N is the last time instant of the conventional STC. It should be noted that the elements of θ̂i(k) in (91) would not be independent with each other, when N-k1+1 is not large enough. Similarly, Qi can be approximated by comparing (88a) with (89a) as follows:(92)Qi≈[Mi(N)εi(N)][Mi(N)εi(N)]T,
where(93)Mi(N)=S̅i(0)ϕi(N)Ri+ϕiT(N)S̅i(0)ϕi(N).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 Qi obtained from (92), while the system is healthy, is so small in general; hence θ̂i(k)≅θ̂i(k-1) in (89a). As a result, it cannot reflect the real parameter variations induced by the unanticipated system faults. Therefore, S̅i(k), Ri, and Qi in (89c) need to be appropriately reset when a fault is detected at time instant kfT. Although the conventional STC with an appropriately reset forgetting factor λ(kf) could improve the estimation of parameter variations, the reset forgetting factor λ(kf) would need some trials for various failure modes. Nevertheless, the resets of S̅i(kf), Qi, and Ri 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 estimateS̅i(k) is given as follows. The variation of the parameter estimate before and after the fault can be approximated as(94)δθ̂i(kf)≈θ̂i(kf)-θ̂i(kh),forkh<kf,
where θ̂i(kh) is the parameter estimate of the healthy system and θ̂i(kf) is the parameter estimate of the faulty system. Then, based on the physical interpretation of (89e), S̅i(k-1) in (89c) can be reasonably reset as(95)S̅i(kf-1)≡diag{diagonal  elements  of[θ̂i(kf)-θ̂i(kh)]×[θ̂i(kf)-θ̂i(kh)]T}≈δ2diag{diagonal  elements  ofθ̂i(kf)θ̂iT(kf)}.
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,Ri in (90) needs to be reset by a moving window with the residual(96)Ri≈1kf-ki+1∑k=kikfεi(k)εiT(k),
where usually 1<(kf-ki)<5. Similarly, Qi in (89c) should be reset by substituting (95) and (96) into (89b) and then (89a) to obtain(97)Qi≈[Mi(kf)εi(kf)][Mi(kf)εi(kf)]T.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 instantkfT of the fault occurrence could be detected by utilizing the ratio of the moving window with the residual Ri in (96) and the average norm of the innovation vectors as (98)RiRfi>γi,
where Rfi=(1/kf)∑k=1kfεi(k)εiT(k) and γi 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 Qi, Ri, and Si(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, resetQi, Ri, and Si(kf) given by (95)–(97).(6)
Go to step(3) and repeat for the modified STC process.


## 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(86)ẋ(t)=f(x)+q(x)u(t)+∑v=1ςβv(t-τv)[fv(x)+qv(u)]+w(t),
where fv(x)+qv(u) 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 βv(t) as shown in Figure 4 [21], where us(t) denotes the unit-step function of time t. fv(x)+qv(u), βv(t), and τv 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) The gradual failure function. (b) The abrupt failure function.
(a)(b)By the optimal linearization method described in Section3.3, the nonlinear faulty system (86) can be accurately linearized as the following linear time-varying model:(87)ẋ(t)=[G(t)+∑v=1ςβv(t-τv)Gv(t)]x(t)+[H(t)+∑v=1ςβv(t-τv)Hv(t)]u(t)+w(t).
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 τv. 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:(88a)θi(k+1)=θi(k)+w̅i(k),(88b)yi(k)=θiT(k)ϕi(k)+v̅i(k),(88c)E̅(w̅i(k)w̅iT(k))=Qi,(88d)E̅(v̅i(k)v̅iT(k))=Ri,where i=1,2,…,p. {w̅i(k)} and {v̅i(k)} are white Gaussian noise sequences with zero means. Applying the Kalman filter to this model (88a) and (88d) gives the following recursive algorithm:(89a)θ̂i(k)=θ̂i(k-1)+Mi(k)εi(k),(89b)Mi(k)=S̅i(k-1)ϕi(k)Ri+ϕiT(k)S̅i(k-1)ϕi(k),(89c)S̅i(k)=S̅i(k-1)-S̅i(k-1)ϕi(k)ϕiT(k)S̅i(k-1)Ri+ϕiT(k)S̅i(k-1)ϕi(k)+Qi,(89d)εi(k)=yi(k)-θ̂iT(k-1)ϕi(k),
with(89e)S̅i(0)=E̅{[θ̂i(0)-θi(0)][θ̂i(0)-θi(0)]T},where S̅i(k) is the covariance matrix of the parameter estimate error. {εi(k)} 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 ofQi, Ri, and S̅i(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 ofQi, Ri, and S̅i(0) to the following moving window-based statistical quantities: (90)Ri≈1N-k1+1∑k=k1Nεi(k)εiT(k),(91)S̅i(0)≈diag{∑k=k1N[(θ̂i(k)-θi(N))(θ̂i(k)-θi(N))T]diagonal  elements  of1N-k1+1×∑k=k1N[(θ̂i(k)-θi(N))(θ̂i(k)-θi(N))T]},
where k1 is the time index after the estimate θ̂i(k) is in steady state, and N is the last time instant of the conventional STC. It should be noted that the elements of θ̂i(k) in (91) would not be independent with each other, when N-k1+1 is not large enough. Similarly, Qi can be approximated by comparing (88a) with (89a) as follows:(92)Qi≈[Mi(N)εi(N)][Mi(N)εi(N)]T,
where(93)Mi(N)=S̅i(0)ϕi(N)Ri+ϕiT(N)S̅i(0)ϕi(N).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 Qi obtained from (92), while the system is healthy, is so small in general; hence θ̂i(k)≅θ̂i(k-1) in (89a). As a result, it cannot reflect the real parameter variations induced by the unanticipated system faults. Therefore, S̅i(k), Ri, and Qi in (89c) need to be appropriately reset when a fault is detected at time instant kfT. Although the conventional STC with an appropriately reset forgetting factor λ(kf) could improve the estimation of parameter variations, the reset forgetting factor λ(kf) would need some trials for various failure modes. Nevertheless, the resets of S̅i(kf), Qi, and Ri 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 estimateS̅i(k) is given as follows. The variation of the parameter estimate before and after the fault can be approximated as(94)δθ̂i(kf)≈θ̂i(kf)-θ̂i(kh),forkh<kf,
where θ̂i(kh) is the parameter estimate of the healthy system and θ̂i(kf) is the parameter estimate of the faulty system. Then, based on the physical interpretation of (89e), S̅i(k-1) in (89c) can be reasonably reset as(95)S̅i(kf-1)≡diag{diagonal  elements  of[θ̂i(kf)-θ̂i(kh)]×[θ̂i(kf)-θ̂i(kh)]T}≈δ2diag{diagonal  elements  ofθ̂i(kf)θ̂iT(kf)}.
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,Ri in (90) needs to be reset by a moving window with the residual(96)Ri≈1kf-ki+1∑k=kikfεi(k)εiT(k),
where usually 1<(kf-ki)<5. Similarly, Qi in (89c) should be reset by substituting (95) and (96) into (89b) and then (89a) to obtain(97)Qi≈[Mi(kf)εi(kf)][Mi(kf)εi(kf)]T.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 instantkfT of the fault occurrence could be detected by utilizing the ratio of the moving window with the residual Ri in (96) and the average norm of the innovation vectors as (98)RiRfi>γi,
where Rfi=(1/kf)∑k=1kfεi(k)εiT(k) and γi 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 Qi, Ri, and Si(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, resetQi, Ri, and Si(kf) 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:(99)y1(k)=y1(k-1)-0.5y1(k-2)+y1(k-2)y2(k-2)+40u1(k-1)-2u2(k-1)+6u12(k-2)+e1(k-1),y2(k)=-0.4y2(k-1)-0.5y1(k-2)-0.2y1(k-2)y2(k-2)cos(y2(k-1))-20u2(k-1)-8u22(k-2)+e2(k-1),
where y1(k) and y2(k) are outputs, u1(k) and u2(k) are inputs, and e1(k) and e2(k) are noises which are zero-mean Gaussian sequences with variances σe12=σe22=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 T is 0.1 second. The identification results of ARMAX and NARMAX models through RELS method are shown in Figures 5 and 6, respectively.The trajectories of the output 1 for ARMAX and NARMAX models.
(a)(b)The trajectories of the output 2 for ARMAX and NARMAX models.
(a)(b)Compare the parameter identification of RELS method based on the intuitive initial parameterθ̂(0)=[I202×9I202×1002×12]T with the OKID-obtained initial parameter θ̂(0) as follows.(1) The Intuitive Initial Parameterθ̂(0) = [I202×9I202×1002×12]T
Set directly the initial parameter of NARMAX model to be[I202×9I202×1002×12]T, 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.The trajectories of outputs using RELS method without OKID.
(a)(b)(c)(d)Figure 8
Convergence of parameter estimation of NARMAX model using RELS method without OKID.(2) The Initial Parameterθ̂(0) Obtained by OKID Method
Assume that the system is excited by the zero mean white noise inputu(t)=[u1(t)u2(t)] with the covariance cov(u)=diag{0.1,0.1}, and the sampling period T is 0.1 second. Then the system realization estimates {Ĝo,Ĥo,Ĉo} and observer gain Ko 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 θ̂(k). 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.The trajectories of outputs using RELS method with OKID.
(a)(b)(c)(d)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 Figure11, where we make the following considerations: (i) the inertial coordinate system {X,Y,Z} is attached to the center of the Earth; (ii) RL∈ℜ3 denotes the position from the origin of {X,Y,Z} to the leader; (iii) the coordinate system {xL,yL,zL} is attached to the leader with the xL axis pointing in the opposite direction as the tangential velocity, the yL axis pointing along the direction of RL, and zL axis being mutually perpendicular to the xLyL-plane; overall, {xL,yL,zL} forms a right-hand coordinate frame; (iv) ρ=[x,y,z]T∈ℜ3 denotes the relative position from the origin of {xL,yL,zL} to the follower.Schematic representations of the MSFF system.
(a)(b)The nonlinear position dynamics of the leader and follower spacecraft with respect to{X,Y,Z} are given by(100)mLR̈L+mL(M+mL)GRL∥RL∥3+FL=uL,mF(R̈L+ρ̈)+mF(M+mF)GRL+ρ∥RL+ρ∥3+FF=uF,
where mL,mF∈ℜ>0 are, respectively, the masses of the leader and follower; FL,FF∈ℜ3 are the disturbance forces, respectively; uL=[uLxuLyuLz]T and uF=[uFxuFyuFz]T are the actual control inputs of the leader and follower, respectively; M=5.974×1024 kg is the mass of the Earth; G=6.673×10-11 N-m2/kg2 is the universal gravity constant.Note that the angular velocity is{0,0,ω}, the relative acceleration is given by ρ̈={ẍ-2ωẏ-ω2x,ÿ+2ωẋ-ω2y,z̈}, and RL={0,∥RL∥,0} is constant in coordinates {xL,yL,zL}. Reflecting M≫mL,mF 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 {xL,yL,zL}, and F=FF-(mF/mL)FL=[FxFyFz]T is the total constant disturbance. Choosing the states variable x=[ρTρ̇T]T=[x1x2x3x4x5x6]T and the control u=uf, the vector field of the concerned system is represented as follows:(101a)ẋ(t)=f(x(t))+B̃u(t)+d(t)=[x4(t)x5(t)x6(t)ω2x1(t)+2ωx5(t)-MGx1(t)∥RL+ρ∥3ω2x2(t)-2ωx4(t)-MG(x2(t)+∥RL∥)∥RL+ρ∥3-MGx3(t)∥RL+ρ∥3]+[031mFI3]u(t)+[000-1mLuLx-1mFFxMG∥RL∥2-1mLuLy-1mFFy-1mLuLz-1mFFz],(101b)y(t)=[y1(t)y2(t)y3(t)]T=Cx(t),C=[I303].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 placedxL-axis thrust denoted by uFxNonideal is illustrated in Figure 12, which results in harmful influence on the other axes thrust actions. Putting the relationship u=uFNonideal into (101a) alters input matrix B=B̃Θ where(102)Θ=[sin(ϕx)cos(θx)sin(ϕy)cos(θy)sin(ϕz)cos(θz)sin(ϕx)sin(θx)sin(ϕy)sin(θy)sin(ϕz)sin(θz)cos(ϕx)cos(ϕy)cos(ϕz)].Figure 12
Nonideal placement ofufx thruster nozzle.Output responses of Scenario1 without fault-tolerant control.
(a)(b)(c)Output responses of Scenario1 with fault-tolerant control.
(a)(b)(c)PWM control inputs of Scenario1 with fault-tolerant control.
(a)(b)(c)Consider the MSFF system with parameters as follows:mL = 1550 kg, mF = 410 kg, RL = [04.224×1070]T m, ω = 7.272×10-5 rad/s, F = 10-5 ×[-1.0256.248-2.415]T N, and uL = 03×1 N; the relative position and velocity are initialized to ρ(0) = [00-200]T m and ρ̇(0) = 03×1 m/s, respectively; the desired relative trajectory is selected as(103)r(t)=[100cos(4ωt)(1-e-0.05t3)100sin(4ωt)(1-e-0.05t3)0].
Besides, we further introduce the nonideal thruster nozzle placements on all axes by assuming ϕx=ϕy=θy=0.5π, θx=0,ϕz=0.05π, and θz=0.1π rad. Let the sampling period T 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, θ̂i(0)i=1,2,3, of the chosen NARMAX model are obtained by OKID method. If a fault is detected at t=kf, matrices Qi, Ri, and S̅i(kf) with the resetting parameter δ=1 in (95) are automatically reset. The fault detection threshold is γ1=γ2=3.5 and γ3=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(104)yi(k)=Fi(Ỹ)+Bi(Ỹ)U+Di(Ỹ)E+εi(k),i=1,2,3,
where(105)Ỹ=[y1𝒜y2𝒜y3𝒜y1ℬy2ℬy3ℬ]T,F1(Ỹ)=a11y1(k-1)+a12y2(k-1)+a13y3(k-1)+a14y1(k-2)+a15y2(k-2)+a16y3(k-2)+a17y12(k-1)+a18y12(k-2)+a19y1(k-1)y2(k-1)+a1,10y1(k-1)y3(k-1)+a1,11y1(k-2)y2(k-1)+a1,12y1(k-2)y3(k-1)+a1,13y1(k-1)y2(k-2)+a1,14y1(k-1)y3(k-2)+a1,15y1(k-2)y2(k-2)+a1,16y1(k-2)y3(k-2)+a1,17y1(k-1)y1(k-2),F2(Ỹ)=a21y1(k-1)+a22y2(k-1)+a23y3(k-1)+a24y1(k-2)+a25y2(k-2)+a26y3(k-2)+a27y22(k-1)+a28y22(k-2)+a29y2(k-1)y1(k-1)+a2,10y2(k-1)y3(k-1)+a2,11y2(k-2)y1(k-1)+a2,12y2(k-2)y3(k-1)+a2,13y2(k-1)y1(k-2)+a2,14y2(k-1)y3(k-2)+a2,15y2(k-2)y1(k-2)+a2,16y2(k-2)y3(k-2)+a2,17y2(k-1)y2(k-2),F3(Ỹ)=a31y1(k-1)+a32y2(k-1)+a33y3(k-1)+a34y1(k-2)+a35y2(k-2)+a36y3(k-2)+a37y32(k-1)+a38y32(k-2)+a39y3(k-1)y1(k-1)+a3,10y3(k-1)y2(k-1)+a3,11y3(k-2)y1(k-1)+a3,12y3(k-2)y2(k-1)+a3,13y3(k-1)y1(k-2)+a3,14y3(k-1)y2(k-2)+a3,15y3(k-2)y1(k-2)+a3,16y3(k-2)y2(k-2)+a3,17y3(k-1)y3(k-2),B1T(Ỹ)=[b11+b17y1(k-1)+b1,10y1(k-2)b12+b18y1(k-1)+b1,11y1(k-2)b13+b19y1(k-1)+b1,12y1(k-2)b14+b1,13y1(k-1)+b1,16y1(k-2)b15+b1,14y1(k-1)+b1,17y1(k-2)b16+b1,15y1(k-1)+b1,18y1(k-2)],B2T(Ỹ)=[b21+b27y2(k-1)+b2,10y2(k-2)b22+b28y2(k-1)+b2,11y2(k-2)b23+b29y2(k-1)+b2,12y2(k-2)b24+b2,13y2(k-1)+b2,16y2(k-2)b25+b2,14y2(k-1)+b2,17y2(k-2)b26+b2,15y2(k-1)+b2,18y2(k-2)],B3T(Ỹ)=[b31+b37y3(k-1)+b3,10y3(k-2)b32+b38y3(k-1)+b3,11y3(k-2)b33+b39y3(k-1)+b3,12y3(k-2)b34+b3,13y3(k-1)+b3,16y3(k-2)b35+b3,14y3(k-1)+b3,17y3(k-2)b36+b3,15y3(k-1)+b3,18y3(k-2)],D1T(Ỹ)=[d11+d17y1(k-1)+d1,10y1(k-2)d12+d18y1(k-1)+d1,11y1(k-2)d13+d19y1(k-1)+d1,12y1(k-2)d14+d1,13y1(k-1)+d1,16y1(k-2)d15+d1,14y1(k-1)+d1,17y1(k-2)d16+d1,15y1(k-1)+d1,18y1(k-2)],D2T(Ỹ)=[d21+d27y2(k-1)+d2,10y2(k-2)d22+d28y2(k-1)+d2,11y2(k-2)d23+d29y2(k-1)+d2,12y2(k-2)d24+d2,13y2(k-1)+d2,16y2(k-2)d25+d2,14y2(k-1)+d2,17y2(k-2)d26+d2,15y2(k-1)+d2,18y2(k-2)],D3T(Ỹ)=[d31+d37y3(k-1)+d3,10y3(k-2)d32+d38y3(k-1)+d3,11y3(k-2)d33+d39y3(k-1)+d3,12y3(k-2)d34+d3,13y3(k-1)+d3,16y3(k-2)d35+d3,14y3(k-1)+d3,17y3(k-2)d36+d3,15y3(k-1)+d3,18y3(k-2)],U=[u1𝒜u2𝒜u3𝒜u1ℬu2ℬu3ℬ]T,E=[ε1𝒜ε2𝒜ε3𝒜ε1ℬε2ℬε3ℬ]T.Simplify the model to form a linear-in-the-parameter expression as follows:(106)yi(k)=θiTϕi(k)+εi(k),i=1,2,3,
where(107)θi=[ai1ai2⋯ai,17bi1bi2⋯bi,18di1di2⋯di,18]T∈ℜ53..ϕi(k) is composed of the corresponding variables associated with each element of the parameter vector θi. ξi=53 is the number of θi. Estimate the parameter θi by using the standard recursive extended-least-squares algorithm, where λ0=λ(0)=0.9,Si(0)=10×Iξi, and the initial parameter θ̂i(0) could be obtained via OKID method.By the optimal linearization method, the optimal linear models ofFi(Ỹ) at sampling time kT are gotten as(108)Fi(Ỹ)=ÃiỸ,i=1,2,3,
where(109)Ãi=[ãi1ãi2ãi3ãi4ãi5ãi6],[ã11ã12ã13ã14ã15ã16]=[a11a12a13a14a15a16]+[2a17a19a1,10a1,17a1,13a1,14a1900a1,1100a1,1000a1,1200a1,17a1,11a1,122a18a1,15a1,16a1,1300a1,1500a1,1400a1,1600]Ỹop-η̃1Ỹop,η̃1=1∥Ỹ∥22×{y1(k-1)×[a17y1(k-1)+a19y2(k-1)+a1,10y3(k-1)+a1,13y2(k-2)+a1,14y3(k-2)+a1,17y1(k-2)]+y1(k-2)×[a18y1(k-2)+a1,11y2(k-1)+a1,12y3(k-1)+a1,15y2(k-2)+a1,16y3(k-2)]}|Ỹ=Ỹop,[ã21ã22ã23ã24ã25ã26]=[a21a22a23a24a25a26]+[0a2900a2,110a292a27a2,10a2,13a2,17a2,140a2,1000a2,1200a2,1300a2,150a2,11a2,17a2,12a2,152a28a2,160a2,1400a2,160]Ỹop-η̃2Ỹop,η̃2=1∥Ỹ∥22×{y2(k-1)×[a29y1(k-1)+a27y2(k-1)+a2,10y3(k-1)+a2,17y2(k-2)+a2,14y3(k-2)+a2,13y1(k-2)]+y2(k-2)×[a2,15y1(k-2)+a2,11y1(k-1)+a2,12y3(k-1)+a28y2(k-2)+a2,16y3(k-2)]}|Ỹ=Ỹop,[ã31ã32ã33ã34ã35ã36]=[a31a32a33a34a35a36]+[00a3900a3,1100a3,1000a3,12a39a3,102a37a3,13a3,14a3,1700a3,1300a3,1500a3,1400a3,16a3,11a3,12a3,17a3,15a3,162a38]Ỹop-η̃3Ỹop,η̃3=1∥Ỹ∥22×{y3(k-1)×[a39y1(k-1)+a3,10y2(k-1)+a37y3(k-1)+a3,14y2(k-2)+a3,17y3(k-2)+a3,13y1(k-2)]+y3(k-2)×[a3,15y1(k-2)+a3,11y1(k-1)+a3,12y2(k-1)+a3,16y2(k-2)+a38y3(k-2)]}|Ỹ=Ỹop.Go(k) and Ho(k) for the associated state-space innovation form are determined to be(110)Go(k)=[-Go1(k)I3-Go2(k)03],Ho(k)=[Ho1(k)Ho2(k)]Co=[I303],
where(111)Go1(k)=[-ã11-ã12-ã13-ã21-ã22-ã23-ã31-ã32-ã33],Go2(k)=[-ã14-ã15-ã16-ã24-ã25-ã26-ã34-ã35-ã36],Ho1(k)=[b11+b17y1𝒜+b1,10y1(k-2)b12+b18y1𝒜+b1,11y1(k-2)b13+b19y1𝒜+b1,12y1(k-2)b21+b27y2𝒜+b2,10y2(k-2)b22+b28y2𝒜+b2,11y2(k-2)b23+b29y2𝒜+b2,12y2(k-2b31+b37y3𝒜+b3,10y3(k-2)b32+b38y3𝒜+b3,11y3(k-2)b33+b39y3𝒜+b3,12y3(k-2)]|Ỹ=Ỹop,Ho2(k)=[b14+b1,13y1𝒜+b1,16y1(k-2)b15+b1,14y1𝒜+b1,17y1(k-2)b16+b1,15y1𝒜+b1,18y1(k-2)b24+b2,13y2𝒜+b2,16y2(k-2)b25+b2,14y2𝒜+b2,17y2(k-2)b26+b2,15y2𝒜+b2,18y2(k-2)b34+b3,13y3𝒜+b3,16y3(k-2)b35+b3,14y3𝒜+b3,17y3(k-2)b36+b3,15y3𝒜+b3,18y3(k-2)]|Ỹ=Ỹop.The weighting matrices are Q=108×diag{7000,7000,6}, R=I3 in (5) and Qo=105×I6,Ro=I3 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,θ̂i(0) for i=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 function2.5+(t/26)(1-2e-6(t-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.Output responses of Scenario2 without fault-tolerant control.
(a)(b)(c)Output responses of Scenario2 with fault-tolerant control.
(a)(b)(c)PWM control inputs of Scenario2 with fault-tolerant control.
(a)(b)(c)


## 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:(99)y1(k)=y1(k-1)-0.5y1(k-2)+y1(k-2)y2(k-2)+40u1(k-1)-2u2(k-1)+6u12(k-2)+e1(k-1),y2(k)=-0.4y2(k-1)-0.5y1(k-2)-0.2y1(k-2)y2(k-2)cos(y2(k-1))-20u2(k-1)-8u22(k-2)+e2(k-1),
where y1(k) and y2(k) are outputs, u1(k) and u2(k) are inputs, and e1(k) and e2(k) are noises which are zero-mean Gaussian sequences with variances σe12=σe22=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 T is 0.1 second. The identification results of ARMAX and NARMAX models through RELS method are shown in Figures 5 and 6, respectively.The trajectories of the output 1 for ARMAX and NARMAX models.
(a)(b)The trajectories of the output 2 for ARMAX and NARMAX models.
(a)(b)Compare the parameter identification of RELS method based on the intuitive initial parameterθ̂(0)=[I202×9I202×1002×12]T with the OKID-obtained initial parameter θ̂(0) as follows.(1) The Intuitive Initial Parameterθ̂(0) = [I202×9I202×1002×12]T
Set directly the initial parameter of NARMAX model to be[I202×9I202×1002×12]T, 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.The trajectories of outputs using RELS method without OKID.
(a)(b)(c)(d)Figure 8
Convergence of parameter estimation of NARMAX model using RELS method without OKID.(2) The Initial Parameterθ̂(0) Obtained by OKID Method
Assume that the system is excited by the zero mean white noise inputu(t)=[u1(t)u2(t)] with the covariance cov(u)=diag{0.1,0.1}, and the sampling period T is 0.1 second. Then the system realization estimates {Ĝo,Ĥo,Ĉo} and observer gain Ko 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 θ̂(k). 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.The trajectories of outputs using RELS method with OKID.
(a)(b)(c)(d)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 Figure11, where we make the following considerations: (i) the inertial coordinate system {X,Y,Z} is attached to the center of the Earth; (ii) RL∈ℜ3 denotes the position from the origin of {X,Y,Z} to the leader; (iii) the coordinate system {xL,yL,zL} is attached to the leader with the xL axis pointing in the opposite direction as the tangential velocity, the yL axis pointing along the direction of RL, and zL axis being mutually perpendicular to the xLyL-plane; overall, {xL,yL,zL} forms a right-hand coordinate frame; (iv) ρ=[x,y,z]T∈ℜ3 denotes the relative position from the origin of {xL,yL,zL} to the follower.Schematic representations of the MSFF system.
(a)(b)The nonlinear position dynamics of the leader and follower spacecraft with respect to{X,Y,Z} are given by(100)mLR̈L+mL(M+mL)GRL∥RL∥3+FL=uL,mF(R̈L+ρ̈)+mF(M+mF)GRL+ρ∥RL+ρ∥3+FF=uF,
where mL,mF∈ℜ>0 are, respectively, the masses of the leader and follower; FL,FF∈ℜ3 are the disturbance forces, respectively; uL=[uLxuLyuLz]T and uF=[uFxuFyuFz]T are the actual control inputs of the leader and follower, respectively; M=5.974×1024 kg is the mass of the Earth; G=6.673×10-11 N-m2/kg2 is the universal gravity constant.Note that the angular velocity is{0,0,ω}, the relative acceleration is given by ρ̈={ẍ-2ωẏ-ω2x,ÿ+2ωẋ-ω2y,z̈}, and RL={0,∥RL∥,0} is constant in coordinates {xL,yL,zL}. Reflecting M≫mL,mF 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 {xL,yL,zL}, and F=FF-(mF/mL)FL=[FxFyFz]T is the total constant disturbance. Choosing the states variable x=[ρTρ̇T]T=[x1x2x3x4x5x6]T and the control u=uf, the vector field of the concerned system is represented as follows:(101a)ẋ(t)=f(x(t))+B̃u(t)+d(t)=[x4(t)x5(t)x6(t)ω2x1(t)+2ωx5(t)-MGx1(t)∥RL+ρ∥3ω2x2(t)-2ωx4(t)-MG(x2(t)+∥RL∥)∥RL+ρ∥3-MGx3(t)∥RL+ρ∥3]+[031mFI3]u(t)+[000-1mLuLx-1mFFxMG∥RL∥2-1mLuLy-1mFFy-1mLuLz-1mFFz],(101b)y(t)=[y1(t)y2(t)y3(t)]T=Cx(t),C=[I303].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 placedxL-axis thrust denoted by uFxNonideal is illustrated in Figure 12, which results in harmful influence on the other axes thrust actions. Putting the relationship u=uFNonideal into (101a) alters input matrix B=B̃Θ where(102)Θ=[sin(ϕx)cos(θx)sin(ϕy)cos(θy)sin(ϕz)cos(θz)sin(ϕx)sin(θx)sin(ϕy)sin(θy)sin(ϕz)sin(θz)cos(ϕx)cos(ϕy)cos(ϕz)].Figure 12
Nonideal placement ofufx thruster nozzle.Output responses of Scenario1 without fault-tolerant control.
(a)(b)(c)Output responses of Scenario1 with fault-tolerant control.
(a)(b)(c)PWM control inputs of Scenario1 with fault-tolerant control.
(a)(b)(c)Consider the MSFF system with parameters as follows:mL = 1550 kg, mF = 410 kg, RL = [04.224×1070]T m, ω = 7.272×10-5 rad/s, F = 10-5 ×[-1.0256.248-2.415]T N, and uL = 03×1 N; the relative position and velocity are initialized to ρ(0) = [00-200]T m and ρ̇(0) = 03×1 m/s, respectively; the desired relative trajectory is selected as(103)r(t)=[100cos(4ωt)(1-e-0.05t3)100sin(4ωt)(1-e-0.05t3)0].
Besides, we further introduce the nonideal thruster nozzle placements on all axes by assuming ϕx=ϕy=θy=0.5π, θx=0,ϕz=0.05π, and θz=0.1π rad. Let the sampling period T 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, θ̂i(0)i=1,2,3, of the chosen NARMAX model are obtained by OKID method. If a fault is detected at t=kf, matrices Qi, Ri, and S̅i(kf) with the resetting parameter δ=1 in (95) are automatically reset. The fault detection threshold is γ1=γ2=3.5 and γ3=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(104)yi(k)=Fi(Ỹ)+Bi(Ỹ)U+Di(Ỹ)E+εi(k),i=1,2,3,
where(105)Ỹ=[y1𝒜y2𝒜y3𝒜y1ℬy2ℬy3ℬ]T,F1(Ỹ)=a11y1(k-1)+a12y2(k-1)+a13y3(k-1)+a14y1(k-2)+a15y2(k-2)+a16y3(k-2)+a17y12(k-1)+a18y12(k-2)+a19y1(k-1)y2(k-1)+a1,10y1(k-1)y3(k-1)+a1,11y1(k-2)y2(k-1)+a1,12y1(k-2)y3(k-1)+a1,13y1(k-1)y2(k-2)+a1,14y1(k-1)y3(k-2)+a1,15y1(k-2)y2(k-2)+a1,16y1(k-2)y3(k-2)+a1,17y1(k-1)y1(k-2),F2(Ỹ)=a21y1(k-1)+a22y2(k-1)+a23y3(k-1)+a24y1(k-2)+a25y2(k-2)+a26y3(k-2)+a27y22(k-1)+a28y22(k-2)+a29y2(k-1)y1(k-1)+a2,10y2(k-1)y3(k-1)+a2,11y2(k-2)y1(k-1)+a2,12y2(k-2)y3(k-1)+a2,13y2(k-1)y1(k-2)+a2,14y2(k-1)y3(k-2)+a2,15y2(k-2)y1(k-2)+a2,16y2(k-2)y3(k-2)+a2,17y2(k-1)y2(k-2),F3(Ỹ)=a31y1(k-1)+a32y2(k-1)+a33y3(k-1)+a34y1(k-2)+a35y2(k-2)+a36y3(k-2)+a37y32(k-1)+a38y32(k-2)+a39y3(k-1)y1(k-1)+a3,10y3(k-1)y2(k-1)+a3,11y3(k-2)y1(k-1)+a3,12y3(k-2)y2(k-1)+a3,13y3(k-1)y1(k-2)+a3,14y3(k-1)y2(k-2)+a3,15y3(k-2)y1(k-2)+a3,16y3(k-2)y2(k-2)+a3,17y3(k-1)y3(k-2),B1T(Ỹ)=[b11+b17y1(k-1)+b1,10y1(k-2)b12+b18y1(k-1)+b1,11y1(k-2)b13+b19y1(k-1)+b1,12y1(k-2)b14+b1,13y1(k-1)+b1,16y1(k-2)b15+b1,14y1(k-1)+b1,17y1(k-2)b16+b1,15y1(k-1)+b1,18y1(k-2)],B2T(Ỹ)=[b21+b27y2(k-1)+b2,10y2(k-2)b22+b28y2(k-1)+b2,11y2(k-2)b23+b29y2(k-1)+b2,12y2(k-2)b24+b2,13y2(k-1)+b2,16y2(k-2)b25+b2,14y2(k-1)+b2,17y2(k-2)b26+b2,15y2(k-1)+b2,18y2(k-2)],B3T(Ỹ)=[b31+b37y3(k-1)+b3,10y3(k-2)b32+b38y3(k-1)+b3,11y3(k-2)b33+b39y3(k-1)+b3,12y3(k-2)b34+b3,13y3(k-1)+b3,16y3(k-2)b35+b3,14y3(k-1)+b3,17y3(k-2)b36+b3,15y3(k-1)+b3,18y3(k-2)],D1T(Ỹ)=[d11+d17y1(k-1)+d1,10y1(k-2)d12+d18y1(k-1)+d1,11y1(k-2)d13+d19y1(k-1)+d1,12y1(k-2)d14+d1,13y1(k-1)+d1,16y1(k-2)d15+d1,14y1(k-1)+d1,17y1(k-2)d16+d1,15y1(k-1)+d1,18y1(k-2)],D2T(Ỹ)=[d21+d27y2(k-1)+d2,10y2(k-2)d22+d28y2(k-1)+d2,11y2(k-2)d23+d29y2(k-1)+d2,12y2(k-2)d24+d2,13y2(k-1)+d2,16y2(k-2)d25+d2,14y2(k-1)+d2,17y2(k-2)d26+d2,15y2(k-1)+d2,18y2(k-2)],D3T(Ỹ)=[d31+d37y3(k-1)+d3,10y3(k-2)d32+d38y3(k-1)+d3,11y3(k-2)d33+d39y3(k-1)+d3,12y3(k-2)d34+d3,13y3(k-1)+d3,16y3(k-2)d35+d3,14y3(k-1)+d3,17y3(k-2)d36+d3,15y3(k-1)+d3,18y3(k-2)],U=[u1𝒜u2𝒜u3𝒜u1ℬu2ℬu3ℬ]T,E=[ε1𝒜ε2𝒜ε3𝒜ε1ℬε2ℬε3ℬ]T.Simplify the model to form a linear-in-the-parameter expression as follows:(106)yi(k)=θiTϕi(k)+εi(k),i=1,2,3,
where(107)θi=[ai1ai2⋯ai,17bi1bi2⋯bi,18di1di2⋯di,18]T∈ℜ53..ϕi(k) is composed of the corresponding variables associated with each element of the parameter vector θi. ξi=53 is the number of θi. Estimate the parameter θi by using the standard recursive extended-least-squares algorithm, where λ0=λ(0)=0.9,Si(0)=10×Iξi, and the initial parameter θ̂i(0) could be obtained via OKID method.By the optimal linearization method, the optimal linear models ofFi(Ỹ) at sampling time kT are gotten as(108)Fi(Ỹ)=ÃiỸ,i=1,2,3,
where(109)Ãi=[ãi1ãi2ãi3ãi4ãi5ãi6],[ã11ã12ã13ã14ã15ã16]=[a11a12a13a14a15a16]+[2a17a19a1,10a1,17a1,13a1,14a1900a1,1100a1,1000a1,1200a1,17a1,11a1,122a18a1,15a1,16a1,1300a1,1500a1,1400a1,1600]Ỹop-η̃1Ỹop,η̃1=1∥Ỹ∥22×{y1(k-1)×[a17y1(k-1)+a19y2(k-1)+a1,10y3(k-1)+a1,13y2(k-2)+a1,14y3(k-2)+a1,17y1(k-2)]+y1(k-2)×[a18y1(k-2)+a1,11y2(k-1)+a1,12y3(k-1)+a1,15y2(k-2)+a1,16y3(k-2)]}|Ỹ=Ỹop,[ã21ã22ã23ã24ã25ã26]=[a21a22a23a24a25a26]+[0a2900a2,110a292a27a2,10a2,13a2,17a2,140a2,1000a2,1200a2,1300a2,150a2,11a2,17a2,12a2,152a28a2,160a2,1400a2,160]Ỹop-η̃2Ỹop,η̃2=1∥Ỹ∥22×{y2(k-1)×[a29y1(k-1)+a27y2(k-1)+a2,10y3(k-1)+a2,17y2(k-2)+a2,14y3(k-2)+a2,13y1(k-2)]+y2(k-2)×[a2,15y1(k-2)+a2,11y1(k-1)+a2,12y3(k-1)+a28y2(k-2)+a2,16y3(k-2)]}|Ỹ=Ỹop,[ã31ã32ã33ã34ã35ã36]=[a31a32a33a34a35a36]+[00a3900a3,1100a3,1000a3,12a39a3,102a37a3,13a3,14a3,1700a3,1300a3,1500a3,1400a3,16a3,11a3,12a3,17a3,15a3,162a38]Ỹop-η̃3Ỹop,η̃3=1∥Ỹ∥22×{y3(k-1)×[a39y1(k-1)+a3,10y2(k-1)+a37y3(k-1)+a3,14y2(k-2)+a3,17y3(k-2)+a3,13y1(k-2)]+y3(k-2)×[a3,15y1(k-2)+a3,11y1(k-1)+a3,12y2(k-1)+a3,16y2(k-2)+a38y3(k-2)]}|Ỹ=Ỹop.Go(k) and Ho(k) for the associated state-space innovation form are determined to be(110)Go(k)=[-Go1(k)I3-Go2(k)03],Ho(k)=[Ho1(k)Ho2(k)]Co=[I303],
where(111)Go1(k)=[-ã11-ã12-ã13-ã21-ã22-ã23-ã31-ã32-ã33],Go2(k)=[-ã14-ã15-ã16-ã24-ã25-ã26-ã34-ã35-ã36],Ho1(k)=[b11+b17y1𝒜+b1,10y1(k-2)b12+b18y1𝒜+b1,11y1(k-2)b13+b19y1𝒜+b1,12y1(k-2)b21+b27y2𝒜+b2,10y2(k-2)b22+b28y2𝒜+b2,11y2(k-2)b23+b29y2𝒜+b2,12y2(k-2b31+b37y3𝒜+b3,10y3(k-2)b32+b38y3𝒜+b3,11y3(k-2)b33+b39y3𝒜+b3,12y3(k-2)]|Ỹ=Ỹop,Ho2(k)=[b14+b1,13y1𝒜+b1,16y1(k-2)b15+b1,14y1𝒜+b1,17y1(k-2)b16+b1,15y1𝒜+b1,18y1(k-2)b24+b2,13y2𝒜+b2,16y2(k-2)b25+b2,14y2𝒜+b2,17y2(k-2)b26+b2,15y2𝒜+b2,18y2(k-2)b34+b3,13y3𝒜+b3,16y3(k-2)b35+b3,14y3𝒜+b3,17y3(k-2)b36+b3,15y3𝒜+b3,18y3(k-2)]|Ỹ=Ỹop.The weighting matrices are Q=108×diag{7000,7000,6}, R=I3 in (5) and Qo=105×I6,Ro=I3 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,θ̂i(0) for i=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 function2.5+(t/26)(1-2e-6(t-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.Output responses of Scenario2 without fault-tolerant control.
(a)(b)(c)Output responses of Scenario2 with fault-tolerant control.
(a)(b)(c)PWM control inputs of Scenario2 with fault-tolerant control.
(a)(b)(c)


## 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θ̂i(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.



---

*Source: 217515-2010-06-28.xml*