An Efficient LQR Design for Discrete-Time Linear Periodic System Based on a Novel Lifting Method

TL;DR

This paper introduces a novel lifting method that transforms discrete-time linear periodic systems into augmented LTI systems, enabling more efficient LQR control design through specialized Riccati equation solutions.

Contribution

A new lifting approach for periodic systems that simplifies LQR design by converting them into augmented LTI systems, improving computational efficiency.

Findings

The proposed method outperforms existing algorithms in computational efficiency.

Simulation results demonstrate successful spacecraft attitude control.

The approach effectively simplifies periodic system control design.

Abstract

This paper proposes a novel lifting method which converts the standard discrete-time linear periodic system to an augmented linear time-invariant system. The linear quadratic optimal control is then based on the solution of the discrete-time algebraic Riccati equation associated with the augmented linear time-invariant model. An efficient algorithm for solving the Riccati equation is derived by using the special structure of the augmented linear time-invariant system. It is shown that the proposed method is very efficient compared to the ones that use algorithms for discrete-time periodic algebraic Riccati equation. The efficiency and effectiveness of the proposed algorithm is demonstrated by the simulation test for the design problem of spacecraft attitude control using magnetic torques.

Tables

Table 1. Table 1: CPU time comparison for problem in [ 16 ]

Samples per period	Algorithm 3.1	Yang Algorithm [16]	Hench-Laub Algorithm [15]
100	0.0097 (s)	0.0757 (s)	0.2711 (s)
500	0.2528 (s)	1.6042 (s)	6.5435 (s)
1000	4.2821 (s)	6.3155 (s)	25.8996 (s)

Table 2. Table 2: CPU time comparison for problem in [ 22 ]

Samples per period	Algorithm 3.1	Yang Algorithm [16]	Hench-Laub Algorithm [15]
100	0.0284 (s)	0.1120 (s)	0.3807 (s)
500	3.6376 (s)	2.5629 (s)	9.0144 (s)
1000	38.4912 (s)	10.0629 (s)	36.0690 (s)

Equations

$${\mathbf{x}}_{k+1}={\mathbf{A}}_k{\mathbf{x}}_k+{\mathbf{B}}_k{\mathbf{u}}_k,$$

$$\underset{N\to \infty }{\lim }\left(\min \frac{1}{2}{\mathbf{x}}_N^{\mathrm{T}}{\mathbf{Q}}_N{\mathbf{x}}_N+\frac{1}{2}\sum _{k=0}^{N-1}{\mathbf{x}}_k^{\mathrm{T}}{\mathbf{Q}}_k{\mathbf{x}}_k+{\mathbf{u}}_k^{\mathrm{T}}{\mathbf{R}}_k{\mathbf{u}}_k\right)$$

$${\mathbf{Q}}_k={\mathbf{Q}}_{k+p}\ge 0,$$

$${\mathbf{R}}_k={\mathbf{R}}_{k+p}>0,$$

$${\mathbf{u}}_k=-({\mathbf{R}}_k+{\mathbf{B}}_k^{\mathrm{T}}{\mathbf{P}}_k{\mathbf{B}}_k{)}^{-1}{\mathbf{B}}_k^{\mathrm{T}}{\mathbf{P}}_k{\mathbf{A}}_k{\mathbf{x}}_k.$$

$${\bf E}_{k}=\left[\begin{array}[]{cc}{\bf I}&{\bf B}_{k}{\bf R}_{k}^{-1}{\bf B}_{k}^{{\rm T}}\\ {\bf 0}&{\bf A}_{k}^{{\rm T}}\end{array}\right]={\bf E}_{k+p},$$

$${\bf F}_{k}=\left[\begin{array}[]{cc}{\bf A}_{k}&{\bf 0}\\ -{\bf Q}_{k}&{\bf I}\end{array}\right]={\bf F}_{k+p}.$$

$${\bf E}_{k}^{-1}=\left[\begin{array}[]{cc}{\bf I}&-{\bf B}_{k}{\bf R}_{k}^{-1}{\bf B}_{k}^{{\rm T}}{\bf A}_{k}^{-{\rm T}}\\ {\bf 0}&{\bf A}^{-{\rm T}}_{k}\end{array}\right]={\bf E}_{k+p}^{-1}.$$

$${\bf F}_{k}^{-1}=\left[\begin{array}[]{cc}{\bf A}_{k}^{-1}&{\bf 0}\\ {\bf Q}_{k}{\bf A}_{k}^{-1}&{\bf I}\end{array}\right]={\bf F}_{k+p}^{-1}.$$

$${\mathbf{\Pi }}_k={\mathbf{E}}_{k+p-1}^{-1}{\mathbf{F}}_{k+p-1}{\mathbf{E}}_{k+p-2}^{-1}{\mathbf{F}}_{k+p-2}\dots {\mathbf{E}}_{k+1}^{-1}{\mathbf{F}}_{k+1}{\mathbf{E}}_k^{-1}{\mathbf{F}}_k={\mathbf{\Pi }}_{k+p},$$

$${\mathbf{\Gamma }}_k={\mathbf{F}}_k^{-1}{\mathbf{E}}_k{\mathbf{F}}_{k+1}^{-1}{\mathbf{E}}_{k+1}\dots ,{\mathbf{F}}_{k+p-2}^{-1}{\mathbf{E}}_{k+p-2}{\mathbf{F}}_{k+p-1}^{-1}{\mathbf{E}}_{k+p-1}={\mathbf{\Gamma }}_{k+p}.$$

$${\mathbf{z}}_{k+p}={\mathbf{\Pi }}_k{\mathbf{z}}_k,$$

$${\mathbf{z}}_k={\mathbf{\Gamma }}_k{\mathbf{z}}_{k+p}.$$

$$\left[\begin{array}[]{cc}{\bf T}_{11k}&{\bf T}_{12k}\\ {\bf T}_{21k}&{\bf T}_{22k}\end{array}\right]^{{\rm T}}\boldsymbol{\Pi}_{k}\left[\begin{array}[]{cc}{\bf T}_{11k}&{\bf T}_{12k}\\ {\bf T}_{21k}&{\bf T}_{22k}\end{array}\right]=\left[\begin{array}[]{cc}{\bf S}_{11k}&{\bf S}_{12k}\\ {\bf 0}&{\bf S}_{22k}\end{array}\right],$$

$${\mathbf{P}}_k={\mathbf{T}}_{21k}{\mathbf{T}}_{11k}^{-1}.$$

$$\left[\begin{array}[]{cc}{\bf W}_{11k}&{\bf W}_{12k}\\ {\bf W}_{21k}&{\bf W}_{22k}\end{array}\right]^{{\rm T}}\boldsymbol{\Gamma}_{k}\left[\begin{array}[]{cc}{\bf W}_{11k}&{\bf W}_{12k}\\ {\bf W}_{21k}&{\bf W}_{22k}\end{array}\right]=\left[\begin{array}[]{cc}{\bf U}_{11k}&{\bf U}_{12k}\\ {\bf 0}&{\bf U}_{22k}\end{array}\right],$$

$${\mathbf{P}}_k={\mathbf{W}}_{21k}{\mathbf{W}}_{11k}^{-1}.$$

$${\mathbf{x}}_1={\mathbf{A}}_0{\mathbf{x}}_0+{\mathbf{B}}_0{\mathbf{u}}_0,$$

$${\mathbf{x}}_2={\mathbf{A}}_1{\mathbf{x}}_1+{\mathbf{B}}_1{\mathbf{u}}_1,$$

$${\mathbf{x}}_3={\mathbf{A}}_2{\mathbf{x}}_2+{\mathbf{B}}_2{\mathbf{u}}_2,$$

$${\mathbf{x}}_4={\mathbf{A}}_0{\mathbf{x}}_3+{\mathbf{B}}_0{\mathbf{u}}_3,$$

$${\mathbf{x}}_5={\mathbf{A}}_1{\mathbf{x}}_4+{\mathbf{B}}_1{\mathbf{u}}_4,$$

$${\mathbf{x}}_6={\mathbf{A}}_2{\mathbf{x}}_5+{\mathbf{B}}_2{\mathbf{u}}_5,$$

$${\mathbf{x}}_7={\mathbf{A}}_0{\mathbf{x}}_6+{\mathbf{B}}_0{\mathbf{u}}_6,$$

$$⋮$$

$${\bar{\mathbf{x}}}_1$$

$${\bar{\mathbf{x}}}_2$$

$${\bar{\mathbf{x}}}_{K+1}$$

$$:=$$

$$\bar{{\bf x}}_{0}=\left[\begin{array}[]{c}{\bf x}_{-2}\\ {\bf x}_{-1}\\ {\bf x}_{0}\end{array}\right]:=\left[\begin{array}[]{c}{\bf 0}\\ {\bf 0}\\ {\bf x}_{0}\end{array}\right],\hskip 7.22743pt\bar{{\bf u}}_{0}=\left[\begin{array}[]{c}{\bf u}_{0}\\ {\bf u}_{1}\\ {\bf u}_{2}\end{array}\right].$$

$$\bar{{\bf x}}_{K}=\left[\begin{array}[]{c}{\bf x}_{p(k-1)+1}\\ {\bf x}_{p(k-1)+2}\\ \vdots\\ {\bf x}_{p(k-1)+p}\end{array}\right],\hskip 7.22743pt\bar{{\bf x}}_{0}:=\left[\begin{array}[]{c}{\bf 0}\\ \vdots\\ {\bf 0}\\ {\bf x}_{0}\end{array}\right],\hskip 7.22743pt\bar{{\bf u}}_{K}=\left[\begin{array}[]{c}{\bf u}_{pk}\\ {\bf u}_{pk+1}\\ \ldots\\ {\bf u}_{pk+p-1}\end{array}\right].$$

$$\bar{{\bf x}}_{K+1}=\left[\begin{array}[]{c}{\bf x}_{pk+1}\\ {\bf x}_{pk+2}\\ \vdots\\ {\bf x}_{pk+p}\end{array}\right],\hskip 7.22743pt\bar{{\bf u}}_{K+1}=\left[\begin{array}[]{c}{\bf u}_{p(k+1)}\\ {\bf u}_{p(k+1)+1}\\ \ldots\\ {\bf u}_{p(k+1)+p-1}\end{array}\right].$$

$${\mathbf{x}}_{pk+1}$$

$${\mathbf{x}}_{pk+2}$$

$${\mathbf{x}}_{pk+p}$$

$$\displaystyle\bar{{\bf x}}_{K+1}:=\left[\begin{array}[]{c}{\bf x}_{pk+1}\\ {\bf x}_{pk+2}\\ \vdots\\ {\bf x}_{pk+p}\end{array}\right]$$

Full text

Keywords: Linear periodic discrete-time system, periodic algebraic Riccati equation, LQR, spacecraft attitude control, magnetic torque.

1 Introduction

Many engineering systems are naturally periodic, for example, spacecraft attitude control using magnetic torques [1], helicopter rotors control system [2], wind turbine control system [3], networked control system [4], and multirate sampled data system [5]. It has been known for about six decades that linear periodic time-varying system can be converted to some equivalent linear time-invariant systems [6, 7]. The most popular and widely used methods that convert the linear periodic time-varying model into linear time-invariant models are the so-called lifting methods proposed in [8, 9]. Although these reduced linear time-invariant models are nice for analysis but they are not very convenient for control system design. For example, the Linear Quadratic Regulator (LQR) design for linear periodic system has been focused on the periodic system not on the equivalent linear time-invariant systems proposed in [8, 9]. This strategy leads to extensive research on the solutions of the periodic Riccati equations (see [10, 11, 12, 13, 14] and references therein). For the discrete-time linear periodic system, two efficient algorithms for Discrete-time Periodic Algebraic Riccati Equation (DPARE) are emerged [15, 16].

In this paper, we propose a novel lifting method that converts the linear periodic system to an augmented Linear Time-Invariant (LTI) system. We show that the LQR design method can be directly applied to this LTI system. Moreover, by making full use of the structure of the augmented LTI system, we can derive a very efficient algorithm. We compare the new algorithm to the ones proposed in [15, 16]. In addition to some simple analysis on the efficiency, we demonstrate the efficiency and effectiveness of the new algorithm by the simulation test for the design problems of spacecraft attitude control using magnetic torques.

The remainder of the paper is organized as follows. Section 2 briefly summarizes the algorithms of [15, 16] so that we can compare the proposed algorithm to the existing ones and analyze the efficiency of these algorithms. Section 3 proposes a novel lifting method and applies some standard discrete-time algebraic Riccati equation result to the augmented LTI model. This leads to a very efficient algorithm for the LQR design for the linear periodic system. Section 4 demonstrates the efficiency and effectiveness of the algorithm by some numerical test. Conclusions are summarized in the last section.

2 Periodic LQR design based on linear periodic system

In this section, we briefly review two efficient algorithms for solving DPARE developed in [15, 16]. This will help us later in the comparison of the proposed method to the existing methods.

Let $p$ be an integer representing the total number of samples in one period in a periodic discrete-time system. We consider the following discrete-time linear periodic system given as follows:

[TABLE]

where ${\bf A}_{k}={\bf A}_{k+p}\in{\bf R}^{n\times n}$ and ${\bf B}_{k}={\bf B}_{k+p}\in{\bf R}^{n\times m}$ are periodic time-varying matrices. For this discrete-time linear periodic system (1), the LQR state feedback control is to find the optimal ${\bf u}_{k}$ to minimize the following quadratic cost function

[TABLE]

where

[TABLE]

and the initial condition ${\bf x}_{0}$ is given. It is well-known that the LQR design for problem (1-2) can be solved by using the periodic solution of the discrete-time periodic algebraic Riccati equation [12]. Two efficient algorithms [15, 16] have been developed to solve $p$ $n$-dimensional matrix Riccati equations to find $p$ positive semidefinite matrices ${\bf P}_{k}$, $k=1,\ldots,p$. Given, ${\bf P}_{k}$, the periodic feedback controllers are given by the following equations:

[TABLE]

We summarize these two algorithms as follows: Let

[TABLE]

[TABLE]

If ${\bf A}_{k}$ is invertible, then ${\bf E}_{k}$ and ${\bf F}_{k}$ are invertible, and

[TABLE]

and

[TABLE]

Let ${\bf y}_{k}$ be the costate of ${\bf x}_{k}$, ${\bf z}_{k}=[{\bf x}_{k}^{{\rm T}},{\bf y}_{k}^{{\rm T}}]^{{\rm T}}$, and

[TABLE]

[TABLE]

The solutions of $p$ discrete-time periodic algebraic Riccati equations are symmetric positive semi-definite matrices, ${\bf P}_{k}$, $k=1,\ldots,p$, which are related to the solutions of either one of the two linear systems of equations [15, 16]:

[TABLE]

[TABLE]

Therefore, ${\bf P}_{k}$, $k=1,\ldots,p$, can be obtained by two methods. The first method uses Schur decomposition:

[TABLE]

where ${\bf S}_{11k}$ is upper-triangular and has all of its eigenvalues inside the unique circle. The periodic solution ${\bf P}_{k}$, $k=1,\ldots,p$, is given by [15]

[TABLE]

The second method uses Schur decomposition:

[TABLE]

where ${\bf W}_{11k}$ is upper-triangular and has all of its eigenvalues outside the unique circle. The periodic solution ${\bf P}_{k}$, $k=1,\ldots,p$, is given by [16]

[TABLE]

Remark 2.1

When ${\bf A}_{k}$ and ${\bf Q}_{k}$ are constant matrices, the second method is much efficient because ${\bf F}_{k}$ becomes a constant matrix and ${\bf F}_{k}^{-1}=\cdots={\bf F}_{k+p-1}^{-1}={\bf F}^{-1}$, which makes the computation of (9) much more efficient than the computation of (8).

3 Periodic LQR design based on linear time-invariant system

We propose a lifting method in this section to convert the discrete-time linear periodic system into an augmented linear time-invariant system. Thereby, the periodic LQR design is reduced to the LQR design for the augmented linear time-invariant system.

To simplify our discussion, let us consider a periodic system with $p=3$. We will use $k$ for the discrete-time in the periodic system and $K$ for the discrete-time in the augmented system.

[TABLE]

We can easily regroup the periodic system and rewrite it as the following form:

[TABLE]

[TABLE]

in general, for $k\geq 0$ ($K\geq 0$), we have

[TABLE]

where

[TABLE]

It is worthwhile to note that (65) is a linear time-invariant system. We can easily extend the result to the general case. Let

[TABLE]

We will use the following notation:

[TABLE]

Theorem 3.1

Given a linear periodic discrete-time system with period of $p$ as follows:

[TABLE]

Then, this discrete-time periodic system is equivalent to the linear time-invariant system given as follows:

[TABLE]

*where $\bar{{\bf A}}\in{\bf R}^{pn\times pn}$ and $\bar{{\bf B}}\in{\bf R}^{pn\times pm}$. Moreover, the structure of $\bar{{\bf B}}$ matrix guarantees the causality of the system (89).       *

It is worthwhile to emphasize that there is no overlap between $\bar{{\bf x}}_{K+1}$ and $\bar{{\bf x}}_{K}$; in addition, there is no overlap between $\bar{{\bf u}}_{K+1}$ and $\bar{{\bf u}}_{K}$. This is the major difference between the proposed lifting method and the existing lifting methods in [8, 9] (see also [14]). This feature makes it possible to apply existing design methods for linear time-invariant system (89) which is equivalent to the linear periodic system (66). In the remainder of the paper, we will investigate the LQR design for the system (89). The LQR state feedback control is to find the optimal $\bar{{\bf u}}_{K}$ to minimize the following quadratic cost function

[TABLE]

where

[TABLE]

are constant matrices and the initial condition $\bar{{\bf x}}_{0}$ is given. It is straightforward to see that the optimal control problem described by (89) and (90) is time-invariant but equivalent to the time-varying periodic system described by (1) and (2). Moreover, the optimal feedback matrix of the system (89 -90) is given as follows:

[TABLE]

where $\bar{{\bf P}}$ is the solution of the following time-invariant algebraic Riccati equation:

[TABLE]

Notice that $\bar{{\bf A}}$ is not invertible, this algebraic Riccati equation cannot be directly solved by using the algorithms for time-invariant algebraic Riccati equation proposed in [17, 18], but it can be solved using the algorithm proposed in [19]. In this paper, we propose a more efficient algorithm than the one of [19]. The new algorithm makes full use of the specific structure of $\bar{{\bf A}}$ in which the first $(p-1)n$ columns are zeros. Denote

[TABLE]

where $\bar{{\bf Q}}_{1}={\rm diag}({{\bf Q}}_{0},\ldots,{{\bf Q}}_{p-2})\in{\bf R}^{(p-1)n\times(p-1)n}$ and $\bar{{\bf Q}}_{2}={{\bf Q}}_{p-1}\in{\bf R}^{n\times n}$,

[TABLE]

where $\bar{{\bf R}}_{1}={\rm diag}({{\bf R}}_{0},\ldots,{{\bf R}}_{p-2})\in{\bf R}^{(p-1)m\times(p-1)m}$ and $\bar{{\bf R}}_{2}={{\bf R}}_{p-1}\in{\bf R}^{m\times m}$. Let

[TABLE]

where $\bar{{\bf A}}_{1}\in{\bf R}^{(p-1)n\times n}$ $\bar{{\bf A}}_{2}\in{\bf R}^{n\times n}$, $\bar{{\bf F}}=\left[\bar{{\bf A}}_{1}^{{\rm T}},\bar{{\bf A}}_{2}^{{\rm T}}\right]^{{\rm T}}\in{\bf R}^{pn\times n}$,

[TABLE]

where $\bar{{\bf B}}_{1}\in{\bf R}^{(p-1)n\times pm}$ and $\bar{{\bf B}}_{2}\in{\bf R}^{n\times pm}$, and

[TABLE]

where $\bar{{\bf P}}_{11}\in{\bf R}^{(p-1)n\times(p-1)n}$, $\bar{{\bf P}}_{12}\in{\bf R}^{(p-1)n\times n}$, $\bar{{\bf P}}_{21}\in{\bf R}^{n\times(p-1)n}$, and $\bar{{\bf P}}_{22}\in{\bf R}^{n\times n}$. Let

[TABLE]

Substituting (94), (95), (103), (114), (115), and (116) into (93) yields

[TABLE]

or equivalently

[TABLE]

This proves $\bar{{\bf P}}_{12}=\bar{{\bf P}}_{21}^{{\rm T}}={\bf 0}$ and $\bar{{\bf P}}_{11}=\bar{{\bf P}}_{11}^{{\rm T}}=\bar{{\bf Q}}_{1}$. By examining the lower right block of (118), we have

[TABLE]

and

[TABLE]

Let

[TABLE]

We can rewrite the lower right block of (118) as follows:

[TABLE]

The Riccati equation (131) is a special case discussed in [20]. An efficient Matlab function dare that implements an algorithm of [20] is available to solve (131).

Remark 3.1

Comparing to the methods described in the previous section which need to solve $p$ $n$-dimensional discrete-time Riccati equations, we need only to solve one $n$-dimensional discrete-time Riccati equation using the method proposed in this section.

To compare the efficiency of the method to the ones developed in [15, 16], we will not use the Matlab function dare because dare calculates more information than the solution of the Riccati equation (131). Let $\tilde{{\bf B}}=\hat{{\bf B}}$, $\tilde{{\bf R}}=\hat{{\bf R}}$,

[TABLE]

and

[TABLE]

Riccati equation (131) can be solved by either eigen-decomposition or Schur decomposition for the following generalized eigenvalue problem [20, page 1748, equation (8)]:

[TABLE]

If $\tilde{{\bf A}}$ is invertible, the problem is reduced to solve the following eigenvalue problem [16, equation (30)]:

[TABLE]

Using Schur decomposition for (137), we have

[TABLE]

where ${\bf S}_{11}$ is upper-triangular and has all of its eigenvalues outside the unique circle. The solution of the discrete algebraic Riccati equation (131) is given by

[TABLE]

We summarize the proposed algorithm as follows:

Algorithm 3.1

Data: ${\bf A}_{0},\ldots,{\bf A}_{p-1}$, ${\bf B}_{0},\ldots,{\bf B}_{p-1}$, ${\bf Q}_{0},\ldots,{\bf Q}_{p-1}$, ${\bf R}_{0},\ldots,{\bf R}_{p-1}$.

• Step 1: Form

[TABLE]

• Step 2: Form $\hat{{\bf A}}$, $\hat{{\bf B}}$, $\hat{{\bf Q}}$, $\hat{{\bf R}}$, and $\hat{{\bf S}}$ using (130).
• Step 3: Find the solution $\hat{{\bf P}}$ of the discrete-time algebraic Riccati equation (131) using the algorithm of **[20]** implemented as dare or using the algorithm of **[16]**, i.e., (138) and (139).
• Step 4: The solution of the discrete-time algebraic Riccati equation (93) is given by

[TABLE]

Given $\bar{{\bf x}}_{K}$, we can calculate the feedback control by (92). Applying this feedback control to (89), we obtain the next state $\bar{{\bf x}}_{K+1}$.

Remark 3.2

The proposed lift method for discrete linear periodic time-varying system can be used for other design methods such as ${\bf H}_{\infty}$ and robust pole assignment design.

4 Implementation and numerical simulation

In this section, we discuss some details of implementation that will reduce some computation time comparing to the directly implementation described in the previous section. We also report the test of the proposed algorithm for some real world problems discussed in [16, 22]. We compare the test result with the ones obtained in [16].

4.1 Implementation consideration

The most expensive calculations in Algorithm 3.1 are the calculation of $\hat{{\bf Q}}$, $\hat{{\bf R}}$ , and $\hat{{\bf S}}$ in Step 2, and the calculation of $\hat{{\bf R}}^{-1}=\tilde{{\bf R}}^{-1}$ in Step 3. It is easy to check (cf. [21]): (1) direct calculation of $\hat{{\bf Q}}$ requires ${\cal O}(2(p-1)^{2}n^{3})+{\cal O}(2(p-1)n^{3})+{\cal O}(n^{2})$ flops, (2) direct calculation of $\hat{{\bf R}}$ requires ${\cal O}(2p(p-1)^{2}n^{2}m)+{\cal O}(2p^{2}(p-1)nm^{2})+{\cal O}(p^{2}m^{2})$ flops, (3) direct calculation of $\hat{{\bf S}}$ requires ${\cal O}(2(p-1)^{2}n^{3})+{\cal O}(2(p-1)n^{3})$ flops, and (4) directly calculation of $\hat{{\bf R}}^{-1}$ requires ${\cal O}(p^{3}m^{3})$ flops. For extremely large $p$, i.e., very long period of the system, the majority of the computation is the computation of $\hat{{\bf R}}$ and $\hat{{\bf R}}^{-1}$.

Let ${\bf Q}_{A}={\bf Q}_{1}^{\frac{1}{2}}\bar{{\bf A}}_{1}\in{\bf R}^{(p-1)n\times n}$ and ${\bf Q}_{B}={\bf Q}_{1}^{\frac{1}{2}}\bar{{\bf B}}_{1}\in{\bf R}^{(p-1)n\times pm}$. We will use using Matlab notation for sub-matrices. Since $\bar{{\bf Q}}_{1}$, $\bar{{\bf Q}}_{2}$, and $\bar{{\bf R}}$ are positive diagonal matrices, we can calculate $\hat{{\bf Q}}_{1}$, $\hat{{\bf S}}_{2}$, and $\hat{{\bf R}}$ in (130) more efficiently as follows:

• for $i=1:(p-1)n$
• $\hskip 21.68121pt{\bf Q}_{A}(i,:)={\bf Q}_{1}^{\frac{1}{2}}(i,i)\bar{{\bf A}}_{1}(i,:)$;
• end
• $\hat{{\bf Q}}={\bf Q}_{A}^{{\rm T}}{\bf Q}_{A}$
• for $i=1:n$
• $\hskip 21.68121pt\hat{{\bf Q}}(i,i)=\hat{{\bf Q}}(i,i)+\bar{{\bf Q}}_{2}(i,i)$;
• end
• for $i=1:(p-1)n$
• $\hskip 21.68121pt{\bf Q}_{B}(i,:)={\bf Q}_{1}^{\frac{1}{2}}(i,i)\bar{{\bf B}}_{1}(i,:)$;
• end
• $\hat{{\bf R}}={\bf Q}_{B}^{{\rm T}}{\bf Q}_{B}$
• for $i=1:pm$
• $\hskip 21.68121pt\hat{{\bf R}}(i,i)=\hat{{\bf R}}(i,i)+\bar{{\bf R}}(i,i)$;
• end
• $\hat{{\bf S}}={\bf Q}_{A}^{{\rm T}}{\bf Q}_{B}$

It is easy to check (cf. [21]) the flops for the following calculations: (1) the calculation of $\hat{{\bf Q}}$ requires ${\cal O}((p-1)n)+{\cal O}((p-1)n^{2})+{\cal O}(2(p-1)^{2}n^{3})+{\cal O}(n)$ flops, (2) the calculation of $\hat{{\bf R}}$ requires ${\cal O}(p(p-1)nm)+{\cal O}(2p^{2}(p-1)nm^{2})+{\cal O}(pm)$, (3), the calculation of $\hat{{\bf S}}$ requires ${\cal O}(2(p-1)pn^{2}m)$ flops, (4) this does not save the computation of $\hat{{\bf R}}^{-1}$.

4.2 Simulation test for the problem of [16]

The first simulation test problem is the spacecraft attitude control design using only magnetic torques discussed in [16]. The number of states of this system is $n=6$. The number of control inputs of this system is $m=3$. The controllability of this problem is established in [23]. In this simulation test, we use the same discrete-time linear periodic model as in [16] with the same parameters, such as the spacecraft inertia matrix, orbital inclination, orbital altitude, weight matrices ${\bf Q}$ and ${\bf R}$, and the same initial conditions.

Using $p=100$, $p=500$, and $p=1000$, we run all three algorithms (one proposed in this paper, one proposed in [16], and one proposed in [15]) for this design and recorded the CPU times for all three algorithms. The result is presented in Table 1.

Clearly, the proposed method is significantly cheaper than the algorithms in [15, 16].

4.3 Simulation test for the problem of [22]

The second simulation test problem is a combined method for the spacecraft attitude and desaturation control design using both reaction wheels and magnetic torques discussed in [22]. The number of states of this system is $n=9$. The number of control inputs of this system is $m=6$. The controllability of problem is guaranteed because three reaction wheels are assumed to be available.

Using the parameters provided in [22], for $p=100$, $p=500$, and $p=1000$, we obtained solutions for the corresponding algebraic Riccati equations and recorded the CPU times for all three algorithms. The result is presented in Table 2.

For this problem, $m=6$ is twice as large as the previous problem, the proposed algorithm is faster than the algorithms developed in [15, 16] when the total number of samples in one period is moderate ($p=100$ samples per period), but when the total number of samples in one period increases (to $p=500$ or $p=1000$ samples per period), the advantage of the proposed algorithm will be lost because the computation of the inverse of $\tilde{{\bf R}}\in{\bf R}^{6000\times 6000}$ is ${\cal O}(p^{3}m^{3})$ which is very expensive.

5 Conclusion

In this paper, we propose a new lifting method to convert the discrete-time linear periodic system to a discrete-time linear time-invariant system. The LQR design for the discrete-time linear periodic system is then reduced to the LQR design for the discrete-time linear time-invariant system. By applying the standard algorithm to the discrete-time algebraic Riccati equation associated with the augmented LTI system and using the special structure of the augmented LTI system, we derived an efficient algorithm for the LQR design for the discrete time linear periodic system. We show that the new algorithm is very efficient compared to the existing ones for the discrete-time periodic algebraic Riccati equation, in particular when the number of samples in one period is moderate. We demonstrated this property by the numerical test for the design problems of spacecraft attitude control using magnetic torques. The proposed lift method for discrete linear periodic time-varying system can be used for other design methods such as ${\bf H}_{\infty}$ and robust pole assignment designs.