State-Constrained Sub-Optimal Tracking Controller for Continuous-Time Linear Time-Invariant (CT-LTI) Systems and Its Application for DC Motor Servo Systems

Abstract

In this paper, we propose an analytic solution of state-constrained optimal tracking control problems for continuous-time linear time-invariant (CT-LTI) systems that are based on model-based prediction, the quadratic penalty function, and the variational approach. Model-based prediction is a concept taken from model-predictive control (MPC) and this is essential to change the direction of calculation for the solution from backward to forward. The quadratic penalty function plays an important role in deriving the analytic solution since it can transform the problem into a form that does not have inequality constraints. For computational convenience, we also propose a sub-optimal controller derived from the steady-state approximation of the analytic solution and show that the proposed controller satisfies the Lyapunov stability. The main advantage of the proposed controller is that it can be implemented in real time with a lower computational load compared to the implicit MPC. Finally, the simulation results for a DC motor servo system are shown and compared with the results of the direct multi-shooting method and the implicit MPC to verify the effectiveness of the proposed controller.

1. Introduction

Recently, interest has been increasing in control systems that require limitations on the state of the target system. For example, optimal trajectory control for industrial robots [1,2], which limits the workspace for co-work between humans and machines, and optimal powertrain control for hybrid vehicle systems [3,4,5], which has limitations on battery capacity, have become more critical to industry. For analytical and computational convenience, the target system for industrial purposes is often linearized, therefore several studies on linear optimal controllers with state constraints were performed. However, these studies imply large computational loads that make it difficult to implement in real time. The following section explains why these computational loads are caused.

1.1. Solutions and Their Approximations of the Optimal Control Problems

Assume that the target system is a continuous-time linear time-invariant (CT-LTI) system as follows:

$$\dot{x}=Ax+Bu,$$

(1)

$$y=C_{y}x,$$

(2)

where $x\left(t\right)$ is the state of the system, $u\left(t\right)$ is the input of the system, and $y\left(t\right)$ is the output of the system. Let the tracking error $e\left(t\right)$ be

$$e=r-y,$$

(3)

where $r\left(t\right)$ is the reference for the output. Then, the linear-quadratic tracking (LQT) problem can be expressed as follows [6,7,8]:

$$\begin{array}{l}\underset{e,u}{\mathrm{minimize} }J=\frac{1}{2}{\displaystyle {\int }_{t_0}^{t_f}\left(e^{T}Qe+u^{T}Ru\right)dt}\\ \mathrm{subject} \mathrm{to} Ax+Bu-\dot{x}=0,\end{array}$$

(4)

where $Q$ and $R$ are positive definite matrices weighting for the tracking error and the input, $J$ is the cost function, $t_0$ is the initial time, and $t_f$ is the final time. Assume that pair $\left(A,B\right)$ is controllable and pair $\left(A,C_y\right)$ is observable. Then, by the Lagrange multiplier method [6,7,8] and the variational approach (Theorem A1), the LQT problem (4) can be transformed into a problem of finding solutions of the Riccati equation with boundary conditions and an auxiliary dynamic equation related to the reference. The Riccati equation and the auxiliary dynamic equation can be derived as follows:

$$\dot{P}=-PA-A^{T}P+PB{R}^{-1}{B}^{T}P-{C_y}^{T}Q{C}_y,$$

(5)

$$\dot{g}=-\left(A^T-PB{R}^{-1}{B}^T\right)g-{C_y}^{T}Qr,$$

(6)

where $P\left(t\right)$ and $g\left(t\right)$ are the solutions to be determined. These solutions are usually calculated in the backward direction in time or by using iterations since $A-B{R}^{-1}{B}^{T}P$, which is the system matrix of the closed-loop system, must have all the negative real parts of the eigenvalues. However, if there are no constraints on the input and the state, it is well known in control engineering that the solutions of the equations can be approximated to their steady-state values, and the solutions become independent of the time. In the case of input constraints alone, the optimal input can be determined by Pontryagin’s minimum principle [6,7,8,9]. In general, the solutions determined by Pontryagin’s minimum principle are closely related with the solutions of the Equations (5) and (6) and have a form simple enough to be implemented on modern microcontrollers. On the other hand, if there are inequality constraints on the state, the problem statement (4) is not valid, therefore it must be redefined as follows:

$$\begin{array}{l}\underset{e,u}{\mathrm{minimize} }J=\frac{1}{2}{\displaystyle {\int }_{t_0}^{t_f}\left(e^{T}Qe+u^{T}Ru\right)dt}\\ \mathrm{subject} \mathrm{to} \{\begin{array}{l}Ax+Bu-\dot{x}=0\\ C_{h}x+w\le 0\end{array} ,\end{array}$$

(7)

where $C_h$ and $w$ are time-invariant parameters of the inequality constraints. The vector inequality $h\le 0$ means that $h_i\le 0 \left(i=1,\cdots ,n\right)$, where $h={\left[\begin{array}{ccc}{h}_1& \cdots & h_n\end{array}\right]}^T$. In this case, the problem related to the direction of the calculation is hard to avoid.

Since the problem (7) is difficult to solve exactly, many kinds of research were proposed to solve this problem approximately. The first group of research consists of numerical approaches including dynamic programming methods [6,10,11] or direct and indirect methods [11,12,13]. In general, these methods discretize the target system and apply numerical methods. Since this procedure does not change the direction of the calculation, most of these methods need backward calculations or iterations for all the time steps. The second group of research comprises model-predictive control (MPC) methods including implicit MPC [14,15,16,17] and explicit MPC [18,19,20]. The main difference from the first group is that time-forward calculation is possible since this method predicts the optimal states and inputs of the target system for a short time ahead. However, this method also requires repetitive calculations for the predictions, and the precision of the calculation decreases if the time length of the prediction is not long enough, therefore it still requires many computations in general. To reduce the computational load, explicit MPC was proposed. For explicit MPC design, methods dividing the state space are essential, but the results of these methods are not easy to analyze in practice [21] since they are based on numerical iterations. In summary, the computational problem of the state-constrained optimal tracking control is that the direction of the exact computation for the solution is a time-backward calculation which leads to increasing computational loads.

1.2. Outline and Scope of the Paper

In this paper, we propose an analytic solution of state-constrained optimal tracking control problems for CT-LTI systems that are based on the model-based prediction, the quadratic penalty method, and the variational approach in Section 2. The model-based prediction is a concept taken from MPC, and this is essential to change the direction of the calculation for the solution from backward to forward. The quadratic penalty method plays an important role in deriving the analytic solution since it can transform the problem into a form that does not have inequality constraints. For computational convenience, we also propose a sub-optimal controller derived from the steady-state approximation of the analytic solution, and show that the proposed controller satisfies the Lyapunov stability in Section 3. Finally, the simulation results for a DC motor servo system are shown and compared with the results of the direct multi-shooting method and implicit MPC to verify the effectiveness of the proposed controller in Section 4.

2. Analytic Solution of State-Constrained Optimal Tracking Problems

In this section, we describe the analytic solution of state-constrained optimal tracking problems. This solution can be derived by using the model-based prediction, inequality constraints using prediction, the quadratic penalty function, and the variational approach.

2.1. Model-Based Prediction

Suppose that the target system is equal to (1), then a model-based prediction with a fixed time interval $\tau$ can be written as follows [16,22]:

$${\widehat{x}}_{\tau }\left(t\right)=x\left(t+\tau |t\right)=e^{A\tau }x\left(t\right)+{\displaystyle {\int }_t^{t+\tau }{e}^{A\left(t+\tau -\eta \right)}Bu\left(\eta \right)d\eta }$$

(8)

Assume that the time interval $\tau$ is short enough to consider the input as a constant. Then, we can approximate (8) as

$${\widehat{x}}_{\tau }\approx A_{d}x+B_{d}u,$$

(9)

where

$$A_d\left(\tau \right)=e^{A\tau } \mathrm{and} B_d\left(\tau \right)={\displaystyle {\int }_0^{\tau }{e}^{A\eta }d\eta \cdot }B.$$

The calculation of the above matrices is described in [22] (pp. 114–117) and this can be performed by using c2d MATLAB^® command, etc.

2.2. Inequality Constraints Using Prediction

Assume that the left sides of the inequality constraints are

$$h_i\left(x\right)=C_{h\_i}x+w_i, \left(i=1,2,\cdots ,n\right),$$

(10)

where

$$C_h=\left[\begin{array}{c}{C}_{h\_1}\\ ⋮\\ C_{h\_n}\end{array}\right] \mathrm{and} w=\left[\begin{array}{c}{w}_1\\ ⋮\\ w_n\end{array}\right].$$

Then, the inequality constraints have a form of

$$h_i\left(x\right)\le 0$$

(11)

Since (11) should be valid at all times, the following inequalities also should be valid:

$$h_i\left({\widehat{x}}_{\tau }\right)\le 0$$

(12)

Let

$$h\left(x,u\right)=\left[\begin{array}{c}{h}_1\left(x\right)\\ ⋮\\ h_n\left(x\right)\\ h_1\left({\widehat{x}}_{\tau }\right)\\ ⋮\\ h_n\left({\widehat{x}}_{\tau }\right)\end{array}\right]=\left[\begin{array}{c}{C}_{h}x+w\\ C_h{A}_{d}x+C_h{B}_{d}u+w\end{array}\right],$$

(13)

where

$$C_h=\left[\begin{array}{c}{C}_{h\_1}\\ C_{h\_2}\\ ⋮\\ C_{h\_n}\end{array}\right], w=\left[\begin{array}{c}{w}_1\\ w_2\\ ⋮\\ w_n\end{array}\right].$$

Then, (11) and (12) can be rewritten by using (13) as follows:

$$h\left(x,u\right)\le 0.$$

(14)

2.3. Quadratic Penalty Function

Suppose that the penalty function $p_i$ of the inequality constraint $h_i\le 0$ is

$$p_i\left(h_i\right)=\frac{1}{2}{\alpha }_i\left(h_i\right)h_i{}^2, \left(i=1,2,\cdots ,2n\right),$$

(15)

where

$${\alpha }_i\left(h_i\right)=\{\begin{array}{l}0, h_i<0\hfill \\ q_i, h_i\ge 0 \hfill \end{array},$$

and $q_i>0$ is the weight for ${\alpha }_i$ [23,24]. As shown in Figure 1, the meaning of the penalty function is the violation costs of the inequality constraint.

The quadratic penalty function is defined as the sum of the penalty functions:

$$p\left(h\right)={\displaystyle \sum _{i=1}^{2n}{p}_i}\left(h_i\right)=\frac{1}{2}{h}^T\left(x,u\right)\mathrm{diag}\left({\alpha }^T\right) h\left(x,u\right),$$

(16)

where

$$\alpha =\left[\begin{array}{c}{\alpha }_x\\ {\alpha }_{\tau }\end{array}\right], {\alpha }_x=\left[\begin{array}{c}{\alpha }_1\\ ⋮\\ {\alpha }_n\end{array}\right], \mathrm{and} {\alpha }_{\tau }=\left[\begin{array}{c}{\alpha }_{n+1}\\ ⋮\\ {\alpha }_{2n}\end{array}\right],$$

and the diag function is defined as

$$\mathrm{diag}\left(\left[\begin{array}{cccc}{\sigma }_1& {\sigma }_2& \cdots & {\sigma }_n\end{array}\right]\right)=\left[\begin{array}{cccc}{\sigma }_1& 0& \cdots & 0\\ 0& {\sigma }_2& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & {\sigma }_n\end{array}\right].$$

2.4. Variational Approach

By using the Lagrange multiplier method and the penalty function method [23,24], Problem (7) can be transformed into the following problem:

$$\underset{x,u}{\mathrm{minimize} }{J}_a={\displaystyle {\int }_{t_0}^{t_f}L\left(x,\dot{x},u\right)dt},$$

(17)

$$L\left(x,\dot{x},u\right)=\frac{1}{2}{\left(r-C_{y}x\right)}^{T}Q\left(r-C_{y}x\right)+\frac{1}{2}{u}^{T}Ru+{\lambda }^T\left(Ax+Bu-\dot{x}\right)+\frac{1}{2}{h}^T\mathrm{diag}\left({\alpha }^T\right) h,$$

(18)

where $\lambda \left(t\right)$ is the Lagrange multiplier. Let the Hamiltonian function of (17) be

$$H\left(x,u\right)=\frac{1}{2}{\left(r-C_{y}x\right)}^{T}Q\left(r-C_{y}x\right)+\frac{1}{2}{u}^{T}Ru+{\lambda }^T\left(Ax+Bu\right)+\frac{1}{2}{h}^T\mathrm{diag}\left({\alpha }^T\right) h,$$

(19)

and $Q_x=\mathrm{diag}\left({\alpha }_x{}^T\right)$, $Q_{\tau }=\mathrm{diag}\left({\alpha }_{\tau }{}^T\right)$. Then, substituting (13) into (19) yields

$$\begin{array}{cc}\hfill H\left(x,u\right)& =\frac{1}{2}{\left(r-C_{y}x\right)}^{T}Q\left(r-C_{y}x\right)+\frac{1}{2}{u}^{T}Ru+{\lambda }^T\left(Ax+Bu\right)+\frac{1}{2}{\left(C_{h}x+w\right)}^T{Q}_x\left(C_{h}x+w\right)\hfill \\ & +\frac{1}{2}{\left(C_h{A}_{d}x+C_h{B}_{d}u+w\right)}^T{Q}_{\tau }\left(C_h{A}_{d}x+C_h{B}_{d}u+w\right).\hfill \end{array}$$

(20)

By Corollary A1, the following equations hold:

$$\begin{array}{cc}\hfill \dot{\lambda }& =C_y{}^{T}Q\left(r-C_{y}x\right)-A^T\lambda -C_h{}^T{Q}_x\left(C_{h}x+w\right)-A_d{}^T{C}_h{}^T{Q}_{\tau }\left(C_h{A}_{d}x+C_h{B}_{d}u+w\right)\hfill \\ & =-A^T\lambda -\left(Q_1+A_d{}^T{C}_h{}^T{Q}_{\tau }{C}_h{A}_d\right)x-A_d{}^T{C}_h{}^T{Q}_{\tau }{C}_h{B}_{d}u+C_y{}^{T}Qr\hfill \\ & -\left(C_h{}^T{Q}_x+A_d{}^T{C}_h{}^T{Q}_{\tau }\right)w,\hfill \end{array}$$

(21)

$$0=u^{T}R+{\lambda }^{T}B+{\left(C_h{A}_{d}x+C_h{B}_{d}u+w\right)}^T{Q}_{\tau }{C}_h{B}_d,$$

(22)

$$Q_1=C_y{}^{T}Q{C}_y+C_h{}^T{Q}_x{C}_h.$$

(23)

Let

$$R_{\tau }=R+B_d{}^T{C}_h{}^T{Q}_{\tau }{C}_h{B}_d.$$

(24)

Then, the optimal input is

$$u=-R_{\tau }{}^{-1}{B}^T\lambda -R_{\tau }{}^{-1}{B}_d{}^T{C}_h{}^T{Q}_{\tau }{C}_h{A}_{d}x-R_{\tau }{}^{-1}{B}_d{}^T{C}_h{}^T{Q}_{\tau }w.$$

(25)

2.5. Analytical Solution of the Problem

The following procedure is the same method used in the derivation of the Riccati Equation (5) and the auxiliary dynamic Equation (6) [7,8].

Theorem 1.

Assume that the costate is proportional to the state. Then, the costate can be written as

$$\lambda =Px-g,$$

(26)

where $P\left(t\right)$ and $g\left(t\right)$are values to be determined. Then, the following dynamic equations hold:

$$\dot{P}=-P{A}_z-A_z{}^{T}P+PB{R}_{\tau }{}^{-1}{B}^{T}P-Q_2,$$

(27)

$$\begin{array}{cc}\hfill \dot{g}& =-\left(A_z{}^T-PB{R}_{\tau }{}^{-1}{B}^T\right)g-C_y{}^{T}Qr\hfill \\ & +\left[C_h{}^T{Q}_x-PB{R}_{\tau }{}^{-1}{B}_d{}^T{C}_h{}^T{Q}_{\tau }+A_d{}^T{C}_h{}^T{Q}_{\tau }\left(I-C_h{B}_d{R}_{\tau }{}^{-1}{B}_d{}^T{C}_h{}^T{Q}_{\tau }\right)\right]w,\hfill \end{array}$$

(28)

where

$$Q_2=Q_1+A_d{}^T{C}_h{}^T{Q}_{\tau }\left(I-C_h{B}_d{R}_{\tau }{}^{-1}{B}_d{}^T{C}_h{}^T{Q}_{\tau }\right)C_h{A}_d,$$

(29)

$$A_z=A-B{R}_{\tau }{}^{-1}{B}_d{}^T{C}_h{}^T{Q}_{\tau }{C}_h{A}_d.$$

(30)

Proof. 

By differentiating (26),

$$\dot{\lambda }=\dot{P}x+P\dot{x}-\dot{g}.$$

(31)

Substituting (1), (25), and (26) into (31) yields

$$\begin{array}{cc}\hfill \dot{\lambda }& =\left(\dot{P}+PA\right)x+PBu-\dot{g}\hfill \\ & =\left(\dot{P}+P{A}_z\right)x-PB{R}_{\tau }{}^{-1}{B}^T\lambda -PB{R}_{\tau }{}^{-1}{B}_d{}^T{C}_h{}^T{Q}_{\tau }w-\dot{g}\hfill \\ & =\left(\dot{P}+P{A}_z-PB{R}_{\tau }{}^{-1}{B}^{T}P\right)x-\dot{g}+PB{R}_{\tau }{}^{-1}{B}^{T}g-PB{R}_{\tau }{}^{-1}{B}_d{}^T{C}_h{}^T{Q}_{\tau }w,\hfill \end{array}$$

(32)

Since the left sides of (21) and (32) are equal, the right sides of (21) and (32) are also equal. Therefore,

$$\begin{array}{cc}\hfill 0& =\left(\dot{P}+P{A}_z-PB{R}_{\tau }{}^{-1}{B}^{T}P+Q_1+A_d{}^T{C}_h{}^T{Q}_{\tau }{C}_h{A}_d\right)x\hfill \\ & -\dot{g}+PB{R}_{\tau }{}^{-1}{B}^{T}g+\left(C_h{}^T{Q}_x+A_d{}^T{C}_h{}^T{Q}_{\tau }-PB{R}_{\tau }{}^{-1}{B}_d{}^T{C}_h{}^T{Q}_{\tau }\right)w-C_y{}^{T}Qr\hfill \\ & +A^T\lambda +A_d{}^T{C}_h{}^T{Q}_{\tau }{C}_h{B}_{d}u.\hfill \end{array}$$

(33)

Substituting (25) and (26) into (33) yields

$$\begin{array}{c}0=\left(\dot{P}+P{A}_z+A_z{}^{T}P-PB{R}_{\tau }{}^{-1}{B}^{T}P+Q_2\right)x-\dot{g}-\left(A_z{}^T-PB{R}_{\tau }{}^{-1}{B}^T\right)g-C_y{}^{T}Qr\\ +\left[C_h{}^T{Q}_x-PB{R}_{\tau }{}^{-1}{B}_d{}^T{C}_h{}^T{Q}_{\tau }+A_d{}^T{C}_h{}^T{Q}_{\tau }\left(I-C_h{B}_d{R}_{\tau }{}^{-1}{B}_d{}^T{C}_h{}^T{Q}_{\tau }\right)\right]w.\end{array}$$

(34)

Since (34) should be valid for all the states, Equations (27) and (28) hold. □

3. State-Constrained Sub-Optimal Tracking Controller

Exact solutions of (27) and (28) have to be calculated in the backward direction in time, and this is not proper for real-time implementations. Therefore, in this section, we propose a sub-optimal controller that is stable and proper for real-time implementations.

3.1. State-Constrained Sub-Optimal Tracking Controller

The steady-state values of (27) and (28) are

$$0=P_s{A}_z+A_z{}^T{P}_s-P_{s}B{R}_{\tau }{}^{-1}{B}^T{P}_s+Q_2,$$

(35)

$$g_s=-{\left(A_c{}^T\right)}^{-1}{C}_y{}^{T}Qr+{\left(A_c{}^T\right)}^{-1}\left[\begin{array}{l}{C}_h{}^T{Q}_x-P_{s}B{R}_{\tau }{}^{-1}{B}_d{}^T{C}_h{}^T{Q}_{\tau }\\ +A_d{}^T{C}_h{}^T{Q}_{\tau }\left(I-C_h{B}_d{R}_{\tau }{}^{-1}{B}_d{}^T{C}_h{}^T{Q}_{\tau }\right)\end{array}\right]w,$$

(36)

where

$$A_c=A_z-B{R}_{\tau }{}^{-1}{B}^T{P}_s.$$

(37)

Substituting (26) into (25) yields

$$u=-\left(R_{\tau }{}^{-1}{B}^T{P}_s+R_{\tau }{}^{-1}{B}_d{}^T{C}_h{}^T{Q}_{\tau }{C}_h{A}_d\right)x+R_{\tau }{}^{-1}{B}^T{g}_s-R_{\tau }{}^{-1}{B}_d{}^T{C}_h{}^T{Q}_{\tau }w.$$

(38)

Notably, (38) becomes the steady-state LQT controller [7,8] that uses steady-state values of (5) and (6) if $Q_x$ and $Q_{\tau }$ are null matrices.

The sub-optimal controller is implemented by calculating (35)–(38), but ${\alpha }_{\tau }$ is needed in these calculations. Since ${\alpha }_x$ indicates current violations of the state constraints and ${\alpha }_{\tau }$ indicates possible violations of the state constraints that occur in a moment, ${\alpha }_{\tau }$ can be identified approximately by the following procedure:

• Identify ${\alpha }_x$ using current state values and calculate $Q_x=\mathrm{diag}\left({\alpha }_x{}^T\right)$.
• Calculate (24) and (35) using an algebraic Riccati equation solver with $Q_{\tau }=0$.
• Calculate (36)–(38) using the result of step 2 and applying $Q_{\tau }=0$.
• Calculate (9) using the result of step 3.
• Identify ${\alpha }_{\tau }$ using the result of step 4.

Since ${\alpha }_x$ has limited values (zero or a fixed value for each element), offline calculation results for lower computational loads can be used in step 2. The sub-optimal controller requires the following procedure in addition to the above procedure.

6.

Calculate $Q_{\tau }=\mathrm{diag}\left({\alpha }_{\tau }{}^T\right)$ using the result of step 5.

7.

Calculate (24) and (35) using an algebraic Riccati equation solver with the result of step 6.

8.

Calculate (36)–(38) using the result of step 7.

Since ${\alpha }_{\tau }$ also has limited values, step 7 can be performed offline. The main advantage of this procedure is that the total computation time of the procedure is limited to a fixed upper boundary. In the case of offline calculations of steps 2 and 7, it is obvious that the total computation load of the proposed controller is lower than that of the implicit MPC.

3.2. Stability of the Proposed Controller

By the matrix inversion lemma [25,26],

$$\begin{array}{cc}\hfill I-C_h{B}_d{R}_{\tau }{}^{-1}{B}_d{}^T{C}_h{}^T{Q}_{\tau }& =I-C_h{B}_d{\left(R+B_d{}^T{C}_h{}^T{Q}_{\tau }{C}_h{B}_d\right)}^{-1}{B}_d{}^T{C}_h{}^T{Q}_{\tau }\hfill \\ & ={\left(I+C_h{B}_d{R}^{-1}{B}_d{}^T{C}_h{}^T{Q}_{\tau }\right)}^{-1},\hfill \end{array}$$

(39)

$$I-Q_{\tau }{C}_h{B}_d{R}_{\tau }{}^{-1}{B}_d{}^T{C}_h{}^T={\left(I+Q_{\tau }{C}_h{B}_d{R}^{-1}{B}_d{}^T{C}_h{}^T\right)}^{-1}.$$

(40)

Therefore, the following equations hold:

$$Q_2={\left[\begin{array}{c}{C}_y\\ C_h\\ C_h{A}_d\end{array}\right]}^T\left[\begin{array}{ccc}Q& 0& 0\\ 0& Q_x& 0\\ 0& 0& Q_{\tau }{\left(I+C_h{B}_d{R}^{-1}{B}_d{}^T{C}_h{}^T{Q}_{\tau }\right)}^{-1}\end{array}\right]\left[\begin{array}{c}{C}_y\\ C_h\\ C_h{A}_d\end{array}\right],$$

(41)

$$Q_{\tau }{\left(I+C_h{B}_d{R}^{-1}{B}_d{}^T{C}_h{}^T{Q}_{\tau }\right)}^{-1}={\left(I+Q_{\tau }{C}_h{B}_d{R}^{-1}{B}_d{}^T{C}_h{}^T\right)}^{-1}{Q}_{\tau }.$$

(42)

These equations show that $Q_2$ is positive semi-definite since $Q$ is positive definite and $Q_x$, $Q_{\tau }$ are positive semi-definite or null. Assume that pair $\left(A_c,\left[\begin{array}{c}{C}_y\\ C_h\\ C_h{A}_d\end{array}\right]\right)$ is observable since $C_y$ is different from $C_h$ in general and suitable $A_d$ may be selected by changing $\tau$. Then, Equation (35) can be rewritten as

$$P_s{A}_c+A_c{}^T{P}_s=-\left(P_{s}B{R}_{\tau }{}^{-1}{B}^T{P}_s+Q_2\right).$$

(43)

Substituting (36)–(38) into (1) and applying $r=0$ and $w=0$ yields

$$\dot{x}=A_{c}x.$$

(44)

Therefore, it is concluded that the closed-loop system is stable by the Lyapunov stability theorem [27,28] if the observability condition is satisfied.

3.3. Model Modification for Input Smoothing

The input generated by the proposed controller may have severe vibration that is not found in numerical solutions. To mitigate this, we propose a modification of the plant model including a low pass filter before the input as shown in Figure 2.

The state equation of the modified plant model is

$$\frac{d}{dt}\left[\begin{array}{c}x\\ u\end{array}\right]=\left[\begin{array}{cc}A& B\\ 0& -\beta I\end{array}\right]\left[\begin{array}{c}x\\ u\end{array}\right]+\left[\begin{array}{c}0\\ \beta I\end{array}\right]\overline{u},$$

(45)

$$y=\left[\begin{array}{cc}{C}_y& 0\end{array}\right]\left[\begin{array}{c}x\\ u\end{array}\right].$$

(46)

Then, the problem (7) should be changed as follows:

$$\begin{array}{l}\underset{e,u,\overline{u}}{\mathrm{minimize} }J=\frac{1}{2}{\displaystyle {\int }_{t_0}^{t_f}\left({\left[\begin{array}{c}e\\ u\end{array}\right]}^T\left[\begin{array}{cc}Q& 0\\ 0& R\end{array}\right]\left[\begin{array}{c}e\\ u\end{array}\right]+{\overline{u}}^T\overline{R}\overline{u}\right)dt}\\ \mathrm{subject} \mathrm{to} \{\begin{array}{l}\left[\begin{array}{cc}A& B\\ 0& -\beta I\end{array}\right]\left[\begin{array}{c}x\\ u\end{array}\right]+\left[\begin{array}{c}0\\ \beta I\end{array}\right]\overline{u}-\left[\begin{array}{c}\dot{x}\\ \dot{u}\end{array}\right]=0\\ \left[\begin{array}{cc}{C}_h& 0\end{array}\right]\left[\begin{array}{c}x\\ u\end{array}\right]+w\le 0\end{array} ,\end{array}$$

(47)

where $\overline{R}$ is the weight for $\overline{u}$. Since this problem can be solved by using the same methods described in Section 2 and Section 3, we omit a detailed description of the solution.

4. Case Study: Application for DC Motor Servo Systems

In this section, to help readers understand how to apply the proposed controller, we show an application of the proposed controller for DC motor servo systems. For precision control of DC motor servo systems, studies including MPC [29], data-driven [30], fuzzy [31], neural network [32], cascade control [33], and the digital twin-based optimization [34] were introduced recently. However, except for MPC, these studies did not consider the state constraints; therefore, we compare MPC and the proposed controller for performance verification. The parameters of the target motor are shown in

Table 1. The parameters of the target motor.

Name	Unit	Value
Rotor inductance ($L_m$)	H	0.0065
Armature resistance ($R_m$)	Ω	2.3
Back EMF constant ($K_b$)	V·sec/rad	0.09
Torque constant ($K_m$)	Nm/A	0.09
Friction coefficient ($B_m$)	Nm·sec	0.00175
Rotor inertia ($J_m$)	Nm·rad	0.0000525
Gear ratio (${\kappa }_m$)		0.25

. The target motor is a 24 V DC brushed gear motor and its rated torque is 2.94 Nm.

Then, the state equation of the target system is

$$\begin{array}{cc}\hfill \dot{x}& =\left[\begin{array}{ccc}-\frac{R_m}{L_m}& -\frac{K_b}{L_m}& 0\\ \frac{K_m}{J_m}& -\frac{B_m}{J_m}& 0\\ 0& {\kappa }_m& 0\end{array}\right]x+\left[\begin{array}{c}\frac{1}{L_m}\\ 0\\ 0\end{array}\right]u\hfill \\ & \approx \left[\begin{array}{ccc}-384.62& -13.85& 0\\ 1714.29& -33.33& 0\\ 0& 0.25& 0\end{array}\right]x+\left[\begin{array}{c}153.85\\ 0\\ 0\end{array}\right]u,\hfill \end{array}$$

(48)

$$y=\left[\begin{array}{ccc}0& 0& 1\end{array}\right]x.$$

(49)

where

$$x=\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right],$$

$x_1$ is the motor current, $x_2$ is the angular speed of the motor, and $x_3$ is the angular position of the motor. Let $\beta =1000$, then the modified plant model is

$$\frac{d}{dt}\left[\begin{array}{c}x\\ u\end{array}\right]=A\left[\begin{array}{c}x\\ u\end{array}\right]+B\overline{u},$$

(50)

$$y=C_y\left[\begin{array}{c}x\\ u\end{array}\right],$$

(51)

where

$$A=\left[\begin{array}{cccc}-384.62& -13.85& 0& 153.85\\ 1714.29& -33.33& 0& 0\\ 0& 0.25& 0& 0\\ 0& 0& 0& -1000\end{array}\right],B=\left[\begin{array}{c}\hfill 0\\ \hfill 0\\ \hfill 0\\ \hfill 1000\end{array}\right],C_y=\left[\begin{array}{cccc}0& 0& 1& 0\end{array}\right]$$

Let $\tau =0.001$, then the parameters for the prediction can be calculated as

$$A_d=\left[\begin{array}{cccc}0.67& -0.01& 0& 0.09\\ 1.41& 0.96& 0& 0.07\\ 0& 0& 1& 0\\ 0& 0& 0& 0.33\end{array}\right],B_d=\left[\begin{array}{c}0.04\\ 0.04\\ 0\hfill \\ 0.67\end{array}\right]$$

The optimal tracking problem is

$$\begin{array}{l}\underset{x,u,\overline{u}}{\mathrm{minimize} }J=\frac{1}{2}{\displaystyle {\int }_{t_0}^{t_f}\left({\left[\begin{array}{c}r-x_3\\ u\end{array}\right]}^T\left[\begin{array}{cc}100& 0\\ 0& 1\end{array}\right]\left[\begin{array}{c}r-x_3\\ u\end{array}\right]+{\overline{u}}^2\right)dt}\\ \mathrm{subject} \mathrm{to} \{\begin{array}{l}A\left[\begin{array}{c}x\\ u\end{array}\right]+B\overline{u}-\left[\begin{array}{c}\dot{x}\\ \dot{u}\end{array}\right]=0\\ C_h\left[\begin{array}{c}x\\ u\end{array}\right]+w\le 0,\end{array}\end{array}$$

(52)

where

$$C_h=\left[\begin{array}{cccc}-1& 0& 0& 0\\ 1& 0& 0& 0\\ 0& -1& 0& 0\\ 0& 1& 0& 0\end{array}\right],w=\left[\begin{array}{c}-3\\ -3\\ -50\\ -50\end{array}\right]$$

In this case, $\left[\begin{array}{c}{C}_y\\ C_h\\ C_h{A}_d\end{array}\right]$ has full rank, therefore the closed-loop system is stable. Available values of $P_s$ in the cases of $Q_{\tau }=0$ are shown in

Table 2. Examples of the solutions of (34).

${\mathsf{\alpha }}_{\mathit{x}}$	${\mathbf{P}}_{\mathit{s}}$	${\mathsf{\alpha }}_{\mathit{x}}$	${\mathbf{P}}_{\mathit{s}}$
$\left[\begin{array}{c}0\\ 0\\ 0\\ 0\end{array}\right]$	$\left[\begin{array}{cccc}0.08& 0.02& 6.50& 0.01\\ *& 0& 1.46& 0\\ *& *& 555.69& 1\\ *& *& *& 0\end{array}\right]$	$\left[\begin{array}{c}10000\\ 0\\ 0\\ 0\end{array}\right]$	$\left[\begin{array}{cccc}12.21& -0.13& 6.63& 1.29\\ *& 0.06& 1.56& -0.04\\ *& *& 561.18& 1\\ *& *& *& 0.20\end{array}\right]$
$\left[\begin{array}{c}0\\ 10000\\ 0\\ 0\end{array}\right]$	$\left[\begin{array}{cccc}12.21& -0.13& 6.63& 1.29\\ *& 0.06& 1.56& -0.04\\ *& *& 561.18& 1\\ *& *& *& 0.20\end{array}\right]$	$\left[\begin{array}{c}0\\ 0\\ 0.001\\ 0\end{array}\right]$	$\left[\begin{array}{cccc}0.08& 0.02& 6.50& 0.01\\ *& 0& 1.46& 0\\ *& *& 555.69& 1\\ *& *& *& 0\end{array}\right]$
$\left[\begin{array}{c}10000\\ 0\\ 0.001\\ 0\end{array}\right]$	$\left[\begin{array}{cccc}12.21& -0.13& 6.63& 1.29\\ *& 0.06& 1.56& -0.04\\ *& *& 561.18& 1\\ *& *& *& 0.20\end{array}\right]$	$\left[\begin{array}{c}0\\ 10000\\ 0.001\\ 0\end{array}\right]$	$\left[\begin{array}{cccc}12.21& -0.13& 6.63& 1.29\\ *& 0.06& 1.56& -0.04\\ *& *& 561.18& 1\\ *& *& *& 0.20\end{array}\right]$
$\left[\begin{array}{c}0\\ 0\\ 0\\ 0.001\end{array}\right]$	$\left[\begin{array}{cccc}0.08& 0.02& 6.50& 0.01\\ *& 0& 1.46& 0\\ *& *& 555.69& 1\\ *& *& *& 0\end{array}\right]$	$\left[\begin{array}{c}10000\\ 0\\ 0\\ 0.001\end{array}\right]$	$\left[\begin{array}{cccc}12.21& -0.13& 6.63& 1.29\\ *& 0.06& 1.56& -0.04\\ *& *& 561.18& 1\\ *& *& *& 0.20\end{array}\right]$
$\left[\begin{array}{c}0\\ 10000\\ 0\\ 0.001\end{array}\right]$	$\left[\begin{array}{cccc}12.21& -0.13& 6.63& 1.29\\ *& 0.06& 1.56& -0.04\\ *& *& 561.18& 1\\ *& *& *& 0.20\end{array}\right]$

. The weight for ${\alpha }_x$ is ${\left[\begin{array}{cccc}10000& 10000& 0.001& 0.001\end{array}\right]}^T$ and $P_s$ was calculated by the MATLAB^® care function. Since the weight for ${\alpha }_{\tau }$ is ${\left[\begin{array}{cccc}100& 100& 10& 10\end{array}\right]}^T$ and similar calculations can be performed in the cases of $Q_{\tau }\ne 0$, the proposed controller can be implemented by these results. Figure 3 shows the simulation results with $r=\pi$ and $x\left(0\right)=0$ without changing $w$. The proposed controller was implemented by using Simulink^® blocks and the ode2 (Heun) fixed-step solver was used for the simulation. The step size was set at 10μs for smooth results. To verify the performance of the proposed controller, we also implemented an implicit MPC by using MATLAB^® Model Predictive Control Toolbox^TM software. A numerical method based on the direct multi-shooting method implemented by using CasADi software [35] was chosen for comparison since it is close to the optimal solution. The sampling time used in both the numerical method and the implicit MPC was set to 100 μs since lower sampling time causes larger computational loads. The prediction and control horizon of the MPC is set to (100, 10).

As shown in Figure 3a–c, the trajectory of the proposed controller and that of the numerical method are very similar, therefore it can be said that the proposed controller is well-approximated. On the other hand, the implicit MPC has different trajectories between 0.2 and 0.4 s though the constraints are maintained. In Figure 3d–e, it seems that the trajectory of the proposed controller has a delay compared to that of the numerical method, but the delay is small enough to be neglected. Since the steady-state LQT controller has similar delay properties [36], it is assumed that the cause of this phenomenon is the steady-state approximation.

Figure 4 shows simulation results with $r=10$, $x\left(0\right)=0$, and $w={\left[\begin{array}{cccc}-10& -10& -100& -100\end{array}\right]}^T$. All other control parameters used in Figure 4 are the same as those used in Figure 3. These results show that the inequality constraints are maintained, and this leads to the conclusion that the proposed controller does not seems sensitive to $r$ or $w$. The results of the implicit MPC in Figure 4a–c also differ from these of other methods therefore it is concluded that the implicit MPC differs from the optimal solution.

Figure 5 shows experimental and simulation results with $r=\pi$, $x\left(0\right)=0$, and $w={\left[\begin{array}{cccc}-3& -3& -50& -50\end{array}\right]}^T$. All other control parameters used in Figure 5 are the same as those used in Figure 3. Texas Instruments LAUNCHXL2-570LC43 and BOOSTXL-DRV8323RS (Dallas, TX, USA) are used to control the target motor, and they are shown in Figure 6b. The target motor is equipped with an incremental encoder as shown in Figure 6a, and the encoder has resolutions of 0.0879 degrees. The controller was designed to have sampling time of 100 μs, implemented by using Simulink^® blocks, and converted to C code by using Embedded Coder^® software and Embedded Coder^® support package for ARM^® Cortex^®-R processors. For comparison, simulations were also performed at the same sampling time in Figure 5. The experimental and simulation results have similar trends in Figure 5, but the experimental results have severe vibrations and performance degradations. Since the real motor has nonlinearities and frictions, it seems that these uncertainties are related to the performance degradations.

5. Discussion

In this paper, we proposed a sub-optimal tracking controller that does not need numerical iterations or backward calculations for state-constrained optimal tracking problems. The main advantage of the proposed controller is that it can be implemented in real time with a lower computational load compared to the implicit MPC. Though there is a delay compared to the results of the numerical method, the simulation results show that the proposed controller has acceptable performance. However, the proposed controller needs to be verified through more application cases including industrial robots, hybrid vehicles, or other control systems. For future work, the proposed controller may be extended to discrete-time systems. Therefore, studies related to the discrete-time optimal tracking controller based on the proposed method are worth researching. In particular, the controller may be applied to path tracking control for autonomous vehicles [37,38], which have recently been the subject of much research.