Application of Terminal Region Enlargement Approach for Discrete Time Quasi Infinite Horizon NMPC

Two norm of a vector $\mathbf{y}$ is defined as $|\mathbf{y}|:=\sqrt{\mathbf{y}^{T}\mathbf{y}}$. Weighted norm having weighting matrix ${\mathbf{A}}$ is defined as $\mathbf{y_{A}}:=\sqrt{\mathbf{y^{T}Ay}}$. $norm(\mathbf{A})$ represents a two norm of matrix $\mathbf{A}$. $\lambda_{min}(\mathbf{A})$ indicates the smallest eigenvalue of $\mathbf{A}$. $\rho_{max}(\mathbf{A})$ indicates the maximum magnitude eigenvalue of $\mathbf{A}$. $\mathbb{R}$ is a set of real numbers. $\mathbb{N}$ is a set of natural numbers.

Consider an autonomous discrete time system

with initial condition ${\bf{x}}(0)={\mathbf{x}}_{0}$. Consider a set $\mathcal{X}\subset\mathbb{R}^{n_{x}}$ such that $\mathbf{0}\in\mathcal{X}$.

For the discrete time system given by (2), QIH-NMPC formulation is stated as follows:

with

subject to

where $\mathbf{W}_{x}$ and $\mathbf{W}_{u}$ are state and input weighting matrices of dimension $\left(n_{x}\times n_{x}\right)$, $\left(n_{u}\times n_{u}\right)$ respectively. $\mathbf{P}$ is the terminal penalty matrix of dimension $\left(n_{x}\times n_{x}\right)$. $\mathbf{W}_{x},\mathbf{W}_{u},\mathbf{P}$ are symmetric positive definite matrices. $N$ is a finite prediction horizon length. $\mathbf{z}(i)$ denotes the predicted state in the QIH-NMPC formulation and $\overline{\mathbf{u}}(i)$ denotes the future input moves. The set $\Omega$ is termed as the terminal region in the neighborhood of the origin. The set $\mathcal{X}_{N}\subset\mathcal{X}\subset\mathbb{R}^{n_{x}}$ is termed as the region of attraction is the set of all feasible initial conditions i.e. it is a set of all initial conditions $\mathbf{x}_{0}$ such that terminal inequality constraint (14) is satisfied with inputs constrained given by equation (9) satisfied.

The terminal region $\Omega$ is chosen as an invariant set for the nonlinear system (2) controlled by local linear controller with gain matrix $\mathbf{K}$. The terminal penalty term is chosen such that for all trajectories starting from any point inside the terminal region $\Omega$, with approximation that a single cost term having larger value that the sum of all the predicted stage cost terms from end of horizon to infinity and is given as follows:

with $\overline{\mathbf{u}}(i)=-\mathbf{K}\mathbf{z}(i)\in\mathcal{U}$ for all $k\geq k+N$ and for all $\mathbf{z}(k+N)\in\Omega$.

It is assumed that the solution to the optimal problem (4) with stage cost defined by (8) with input set given by (9) subject to the predicted state dynamics (10) with initial condition (11), predicted states constraint (12), future input moves constraint (13) and terminal state constraint (14) i.e. $\overline{\mathbf{u}}_{[k,k+N-1]}^{*}$ exists and can be computed numerically. Controller is implemented as a moving horizon framework. Accordingly, only the first control move $\mathbf{u}(k)=\overline{\mathbf{u}}^{*}(k)$ is implemented in the plant. Entire process is repeated at the next sampling instant $k+1$. The term Quasi Infinite is because of the fact that the NMPC formulation depicts the stability properties of the infinite horizon formulation, however, the actual implementation contains finite horizon. Such implementation is achieved with the help of the equation (15). However, only ensuring the terminal penalty term satisfying the condition (15) is not sufficient to guarantee the nominal asymptotic stability of the NMPC controller, hence terminal constraint as given by (14) becomes inevitable. It may be noted that local linear controller with gain matrix $\mathbf{K}$ is not used for implementation of the NMPC controller and is only a mathematical construct to characterize the terminal region $\Omega$.

Before proceeding to application of the proposed arbitrary controller based approach and LQR based approach, a look at Yu et al’s approach from [18] is necessary. Consider, Jacobian linearization of the nonlinear system (2) in the neighborhood the origin as,

where

One additional assumption is required at this stage.

D7

The linearized system (16) is stabilizable.

Yu et al. characterize the terminal region as,

where linear gain $\mathbf{L}$ and the terminal penalty matrix $\mathbf{P}$ are the steady state solutions of the modified Lyapunov equation given as follows:

where $\mathbf{\Phi}_{L}=\mathbf{\Phi-\Gamma L}$ and parameter $\kappa$ is chosen such that $1<\kappa<1/\left[\rho_{max}\left(\mathbf{\Phi}_{L}\right)\right]$ with $\rho_{max}(\mathbf{\Phi}_{L})$ being the maximum amplitude eigenvalue of $\mathbf{\Phi}_{L}$. Since the gain $\mathbf{L}$ is stabilizing, all the eigenvalues of $\mathbf{\Phi}_{L}$ lie inside the unit circle, hence $\rho_{max}(\mathbf{\Phi}_{L})<1$. It can be noted that once stage cost weighting matrices $\mathbf{W}_{x},\mathbf{W}_{u}$ are chosen, there is barely any degree of freedom left to the designer for shaping of the terminal region. This results in very conservative terminal regions. In addition, for the systems where there is uncontrollable eigenvalue of $\mathbf{\Phi}$ is close to unit circle, value of $\kappa$ becomes nearly $1$, resulting in a negligible margin for shaping of the terminal region. These limitations are overcome by using the arbitrary controller based approach wherein additive tuning matrices are introduced for providing large degrees of freedom for enlarging of the terminal region and is presented in the subsequent section.

According to the design of the arbitrary controller based approach, the gain matrix $\mathbf{K}$ can be any arbitrary stabilizing linear controller. However, in order to simplify the computations, simulation results are presented with the following choice. Identical linear controller gain is also used in the case of LQR based approach. Controller gain $\mathbf{K}$ is the steady state solution of the simultaneous equations (50) and (51).

The tuning matrix $\Delta\mathbf{Q}$ for the arbitrary controller based approach can be any positive definite matrix. However, in order to simplify and structure the computations of the terminal region, following parameterization is carried out:

In order to further simplify the numerical computation of the terminal region using arbitrary controller based approach, additional parameterization is carried out and is given as follows:

Note that it is sufficient to have $\widetilde{\mathbf{W}}_{x}>\mathbf{W}_{x}$ or $\widetilde{\mathbf{W}}_{u}>\mathbf{W}_{u}$ to satisfy $\Delta\mathbf{Q}>0$, however, usually both $\widetilde{\mathbf{W}}_{x}>\mathbf{W}_{x}$ and $\widetilde{\mathbf{W}}_{u}>\mathbf{W}_{u}$ is preferred in practice. Using the matrices (53) into (52),

where $\rho_{x}$ and $\rho_{u}$ are the tuning scalars. For the arbitrary controller based approach, conditions $\rho_{x}>0$, and $\rho_{u}>0$ and for the LQR based approach, conditions $\rho_{x}>1$ and $\rho_{u}>1$ are imposed.

Efficacy of having two tuning parameters is efficiently demonstrated using the case study in the next sub-section. In the case study, in the first iteration, parameter $\rho_{x}$ is increased from $0$ (for arbitrary controller based approach) or $1$ (for LQR based approach) to a large value, keeping $\rho_{u}$ constant at $0$ or $1$ respectively. A value of $\rho_{x}$ is chosen which leads to larger terminal region in the first iteration. Subsequently, in the second iteration, parameter $\rho_{u}$ is increased from $0$ (for arbitrary controller based approach) or $1$ (for LQR based approach) to larger value, keeping $\rho_{x}$ as constant to the values obtained in the first iteration.

Discrete time system dynamics are given by (55)-(56). These equations are obtained by using the continuous time equations given in [19, 9] with a simple first order Euler integration to obtain the discrete time equations identical to the ones given in [18].

where parameter $\mu_{0}=0.5$, $\mathbf{x}(k)=\left[x_{1}(k)~{}~{}x_{2}(k)\right]^{T}$ is the state vector, and $u(k)$ is the input. The chosen equilibrium operating point is given as $\mathbf{x}(k)=\left[0~{}~{}0\right]^{T}$ and $\mathbf{u}(k)=\left[0\right]$. It may be noted that sampling interval chosen for computations is $T=0.1$s, which is identical to the one chosen in [18]. The eigenvalues of the linearized matrix at the origin are $(0.9048,1.1052)$, indicating that the chosen operating point is unstable resulting in a challenging control problem.

The discrete time QIH-NMPC weighting matrices are chosen identical to the one given in [18] and are given as follows:

The input constraint used in the QIH-NMPC formulation is identical to the one given in [9, 18] and is given as follows:

Table 1 presents the terminal regions obtained using the arbitrary controller based approach when only one tuning parameter $\rho_{x}$ is varied for the two state system.

From the table 1, it may be observed that the terminal region area has nearly saturated at $\rho_{x}=100$ and is used as the value for the second iteration. Table 2 presents the terminal regions obtained using the arbitrary controller based approach when the second tuning parameter $\rho_{u}$ is varied, keeping first tuning parameter $\rho_{x}$ as constant. It may be noted that in this particular system, value of $\alpha$ which defines the terminal region boundary are found identical to its upper bound $\gamma$ and it probably because of the mild nonlinearity of the system dynamics near the origin.

It can be noticed from the table 2 that the terminal region area has reached at $\rho_{u}=100$.

Similar to the arbitrary controller based approach, LQR based approach is implemented. Table 3 presents the terminal regions obtained using the LQR based approach when only one tuning parameter $\rho_{x}$ is varied and using the inequality method for the two state system.

From the table 3, terminal region with largest area is obtained at $\rho_{x}=50$, which is the value used for the second iteration. Table 4 presents the terminal regions obtained using the LQR based approach when the second tuning parameter $\rho_{u}$ is varied, keeping first tuning parameter $\rho_{x}$ as constant.

It can be noticed from the table 4 that the terminal region area has reached maximum saturation around $\rho_{u}=100$.

Table 5 presents a comparison of the terminal regions obtained using approach given in [18], Arbitrary Controller (AC) based approach, and LQR based approach for the two state system system. Table presents the tuning parameters, values of $\gamma$ (upper bound on $\alpha$) and $\alpha$ (which defines the terminal region boundary) and the area of each of the terminal regions.

From the table 5, it can be observed that the largest terminal regions obtained using the arbitrary controller based approach and LQR based approach are approximately 10.4723 and 9.5055 times larger by area respectively when compared to the largest terminal region obtained using the approach given by Yu et al. in [18]. It also illustrates the significant impact of the second degree of freedom in the form of tuning parameter $\rho_{u}$ (or tuning parameter matrix $\widetilde{\mathbf{W}}_{u}$) on the size of the terminal region.

The terminal region obtained using the arbitrary controller based approach (with $\rho_{x}=100$ and $\rho_{u}=100$) having largest area is $\Omega=\left\{\mathbf{x}^{T}\mathbf{P}\mathbf{x}\leq 106.5\right\}$ with an area of $0.2838$. Corresponding values of linear gain matrix $\mathbf{K}$ and terminal penalty matrix $\mathbf{P}$ are given as follows:

The terminal region obtained using the LQR based approach (with $\rho_{x}=50$ and $\rho_{u}=100$) having largest area is $\Omega=\left\{\mathbf{x}^{T}\mathbf{P}\mathbf{x}\leq 95.677\right\}$ with an area of $0.2576$. Corresponding values of linear gain matrix $\mathbf{K}$ and terminal penalty matrix $\mathbf{P}$ are given as follows:

Note that the gain matrix is identical in both the approaches by choice. Discrete time QIH-NMPC controller simulations are carried out using MATLAB software. Sampling interval used for computer simulations is $0.1$ seconds. Initial condition for the state is identical to the one given in [18] and is given as $\mathbf{x}_{0}=[-3~{}~{}2]^{T}$.

Table 6 gives minimum prediction horizon length required to satisfy the terminal constraint. It can be observed that the minimum prediction horizon required in the case of arbitrary controller based approach and LQR based approach is significantly smaller as compared to the prediction horizon required in the case of terminal constraint obtained using Yu et al.’s approach. It is well established in the literature that the computation time required for NMPC iterations increases exponentially with the prediction horizon length [4, 9, 6]. Results obtained effectively validate the advantage of having larger terminal regions.

Figure 1 depicts the plot of the boundaries of the terminal regions obtained using three approaches namely a) approach given by Yu et al. in [18], b) LQR based approach, and c) Arbitrary controller based approach for the two state system for an identical sampling interval of $0.1$s. It may be noted that the terminal regions obtained using arbitrary controller based approach and LQR based approach are significantly larger than the approach by Yu et al. In addition there are several states which lie inside or near to the larger terminal regions obtained using the proposed approaches resulting in feasible initial conditions with a very small prediction horizon length, which will require a relatively larger prediction horizon lengths when the terminal region given in [18] is used as terminal inequality constraint.

Approaches presented in the literature for the terminal region characterization for the discrete time QIH-NMPC formulations provide limited degrees of freedom and often result in a conservative terminal region. Two approaches are presented in this work namely an arbitrary controller based approach and an LQR based approach. Proposed approaches provide two additive matrices as tuning parameters for enlargement of the terminal region. Both the approaches are scalable to system of any state and input dimension.

Efficacy of the approaches is demonstrated using benchmark two state system. It is observed that the terminal region obtained using the arbitrary controller based approach and LQR based approach are approximately 10.4723 and 9.5055 times larger by area measure when compared to the largest terminal region obtained using the approach given by Yu et al. [18]. Future research would involve choosing a completely arbitrary linear controller gain and parameterizing the approaches using different combinations of additive tuning matrices.