Structured linear quadratic control computations over 2D grids

In this paper, we use $\mathbb{R}:=(-\infty,+\infty)$ and $\mathbb{N}:=\{1,2,\ldots\}$ to denote the set of real numbers and natural numbers, respectively. We denote an $n$-dimensional real vector with $\mathbb{R}^{n}$ and a real matrix with $n$ rows and $m$ columns with $\mathbb{R}^{n\times m}$. $\mathcal{A}\backslash\{b\}$ denotes all elements of the set $\mathcal{A}$ except element $b$. $A^{\prime}$ denotes the transpose of a matrix. $A\succ 0$ means the symmetric matrix $A=A^{\prime}\in\mathbb{R}^{n\times n}$ is positive definite (i.e., there exists $c>0$ such that $x^{\prime}Ax\geq cx^{\prime}x$ for all $x\in\mathbb{R}^{n}$) and $A\succeq 0$ means $A$ is a positive semi-definite (i.e., $x^{\prime}Ax\geq 0$ for all $x\in\mathbb{R}^{n}$). “$\mathrm{blkdiag}(.)$”” represents construction of block diagonal matrix from input arguments and “$\mathrm{col}(.)$” represents concatenation of input arguments as column vector.

We are interested in solving finite-horizon linear-quadratic (LQ) optimal control problems with discrete-time dynamics arising from the interconnection of sub-systems in a two-dimensional (2D) grid graph structure as shown in Fig.1.

Each node is a sub-system and is labeled with subscript $(i,j)$ with $i\in\mathcal{K}:=\{1,2,\ldots,K\}\subset\mathbb{N}$ and $j~{}\in~{}\mathcal{N}:=\{1,2,\ldots,N\}\subset\mathbb{N}$, where $K$ is the total number of sub-systems along the vertical direction and $N$ is the total number of sub-systems in the horizontal direction. We assume that there is a bi-directional flow of information between adjacent sub-systems, as shown in Fig.1. That is, each sub-system $(i,j)$ is coupled, at max, with $4$ of the neighbouring sub-systems. Spatio-temporal structure of this kind is relevant in many applications such as automated water irrigation networks [1], radial power distribution networks [2], control of traffic signals [3], spreading processes such as wildfires [4, 5] and in the discretization of 2D partial differential equations [6].

The LQ optimal control problem that we are interested in is a large-scale quadratic program with O(KN T) decision variables, where T is the length of the time horizon. We can solve this problem by solving the linear system of equations that arise from the corresponding first-order optimality conditions, known as KKT conditions. Due to the special spatio-temporal structure of the problem, this linear system of equations is also highly structured.

Previous research has studied the temporal structure of the LQ optimal control problem in [7, 8, 9, 10], but these approaches do not utilize the structure in the spatial dimensions of the problem. Structured solvers for LQ optimal control problems related to path-graph networks (i.e., 1D chains) are proposed in [11, 12, 13], and the computations are decomposable and amenable to distributed implementation on parallel processors with path-graph data exchange. However, these approaches do not exploit the structure in both spatial dimensions simultaneously, as we do in this paper. In [14, 15], structured methods are proposed for solving optimal control problems of multi-agent networks, which can be used to solve the LQ optimal control problem we are interested in. However, these approaches often take a large number of iterations to converge. All of these approaches offer decomposable computations, but they do not exploit the complete three-dimensional spatio-temporal structure of the problem. In [16], a structured solver is proposed for a special case where sub-systems are connected in a directed tree graph structure.

The main contribution of this paper extends the developments presented in [13] to the case of 2D grids. We present a structured solver based on preconditioned conjugate gradient method (PCGM) for solving the KKT conditions. Specifically, we propose a structured preconditioner based on nested block Jacobi method (NBJM) to improve the conditioning of the linear system. The proposed approach retains structure in all three dimensions of the problem, extending the work in [13] which explored structure in two dimensions. We show that the proposed preconditioner is positive-definite. Moreover, the arithmetic complexity of each PCGM step scales as $O(KNT)$. The computations are decomposeable and amenable to distributed implementation on $O(NT)$ parallel processors connected in a 2D grid with localized data exchange.

The rest of the paper is organized as follows. In Section II, we present the structured formulation of the problem and the corresponding KKT conditions. In Section III, we provide an overview of the PCGM. The structured preconditioner based on NBJM is developed in Section IV. We explore performance of the proposed PCGM numerically, for a 2D mass-spring-damper network and an automated water irrigation network, and present results in Section V. Finally, some concluding remarks and possible directions of future work are presented in VI.

The state space dynamics of each sub-system $(i,j)$ are

where $x_{i,j,t}\in\mathbb{R}^{n_{i,j}}$ and $u_{i,j,t}\in\mathbb{R}^{m_{i,j}}$ are the state and input of sub-system $(i,j)\in\mathcal{K}\times\mathcal{N}$ at time $t\in\mathcal{T}:=\{0,1,\ldots,T\}\subset\mathbb{N}\cup\{0\}$, respectively. Initial conditions are given by $x_{i,j,0}=\gamma_{i,j}\in\mathbb{R}^{n_{i,j}}$ and the spatial boundary conditions are given by $x_{0,j,t}=\underline{\alpha}_{j,t}\in\mathbb{R}^{n_{0,j}}$, $x_{K\!+\!1,j,t}=\overline{\alpha}_{j,t}\in\mathbb{R}^{n_{K\!+\!1},j}$ and $x_{i,0,t}=\underline{\beta}_{i,t}\in\mathbb{R}^{n_{i,0}}$, $x_{i,N\!+\!1,t}=\overline{\beta}_{i,t}\in\mathbb{R}^{n_{i,N\!+\!1}}$ for $t\in\mathcal{T}$.

We are interested in the following finite-horizon linear-quadratic (LQ) optimal control problem:

where $\ell_{i,j,t}(x,u)=x^{\prime}Q_{i,j,t}x+u^{\prime}R_{i,j,t}u$. For $i\in\mathcal{K}$, $j\in\mathcal{N}$ and $t\in\mathcal{T}\backslash\{T\}$, it is assumed that $Q_{i,j,t}=Q_{i,j,t}^{\prime}\succ 0$, and $R_{i,j,t}=R_{i,j,t}^{\prime}\succ 0$. Moreover, for every $i\in\mathcal{K}$, $j\in\mathcal{N}$, $Q_{i,j,T}\succ 0$, but $R_{i,j,T}=0$, so that $u_{i,j,T}$ can be removed as a decision variable.

Note that the problem (2) is a large-scale LQ optimal control problem with $O(KNT)$ decision variables. While the cost is separable across the sub-systems in both dimensions as well as the time horizon, the problem is not separable due the spatial coupling between the states of neighbouring sub-systems, and the inter-temporal coupling in the equality constraint (2b). Solution of problem (2) can obtained by solving the large-scale linear system of equations arising from the corresponding first-order optimality conditions known as KKT conditions. However, due to the structured coupling of variables (i.e., states) in the equality constraint (2b), this linear system of equations is also structured.

The positive definite linear system of equations (7a) can be solved using conjugate gradient method (CGM) [18]. Let $e^{(k)}=\tilde{\delta}^{(k)}-\tilde{\delta}^{*}$ be the error between $k$-th iterate $\tilde{\delta}^{(k)}$ of the CGM and the exact solution $\tilde{\delta}^{*}$ of (7a). It can be shown that $e^{(k)}$ satisfies the following [19, Thm. 6.29]:

where $\|e\|_{\Delta}=e^{\prime}\Delta e$, $\kappa(\Delta)=\lambda_{\mathrm{max}}(\Delta)\big{/}\lambda_{\mathrm{min}}(\Delta)$ is the condition number and $\lambda_{\mathrm{max}}(\Delta)$, $\lambda_{\mathrm{min}}(\Delta)$) are the maximum and minimum eigenvalue of $\Delta$, respectively. It shows that the CGM converges faster if $\kappa(\Delta)$ is closer to $1$. Therefore, to improve the performance of CGM, the system of equations (7a) is transformed such that the condition number of the transformed coefficient matrix is improved. The transformed system of equations is known as preconditioned system. Let $P=P^{\prime}$ be a positive definite matrix, then solving (7a) is equivalent to solving the transformed system

where $\tilde{\gamma}=P^{1/2}\tilde{\delta}$ and $r_{\tilde{\gamma}}=P^{-1/2}r_{\tilde{\delta}}$. An efficient implementation of preconditioned CGM (PCGM) is given in Algorithm 1 from [20].

Note that lines 7 to 14 constitute one PCG step. The transformed system $P^{-1/2}\Delta P^{-1/2}$ is never formed explicitly. Instead, the preconditioning is performed by solving a linear system of equations at lines 3 and 12 of Algorithm 1. Choice of preconditioner $P$ determines the overall complexity of PCGM. Ideally, if we use the preconditioner $P=\Delta$, this will result in $P^{-1/2}\Delta P^{-1/2}=I$. However, the linear system to be solved at the preconditioning lines 3 and 12 becomes the original problem. The incomplete sparse LU factorization of the coefficient matrix $\Delta$ can be used as a preconditioner [21, 22]. However, to ensure that the resulting preconditioner is both positive definite and effective, it may be necessary to use incomplete LU factors that are denser (i.e., have less structure) than $\Delta$.

In the next section, we devise a structured preconditioning approach based on nested block Jacobi iterations. We apply a fixed number of inner iterations at each outer iteration of block Jacobi method to approximately implement lines 3 and 12 with $P=\Delta$. The proposed approach is inspired by the ideas from [23, 24, 25, 26, 27, 28]. It extends the developments presented in [13] where same ideas have already been explored within the context of 1D chains.

Let $\Delta=\tilde{\Psi}-\tilde{\Xi}\succ 0$ is a matrix splitting such that

Then, the block Jacobi method for solving (7a) involves the following:

where $s$ is the number of block Jacobi iterations. These iterations converge to the solution of (7a) if and only if

where $\rho(\cdot)$ denotes spectral radius [20, Thm. 2.16]. Given that $\Delta$ is a positive definite block tri-diagonal matrix, and the splitting (14), condition (16) always holds when (15) is solved exactly [20, Lem. 4.7 with Thm. 4.18]. However, instead of solving (15) exactly, we propose a fixed number of inner iterations. It is apparent from (11), that $\tilde{\Psi}\succ 0$ has a block tri-diagonal structure with $\hat{\Omega}_{1,t}=0$ and $\hat{\Omega}_{V+1,t}=0$ for $t\in\mathcal{T}$. Let $\tilde{\Psi}=\tilde{\Phi}-\tilde{\Omega}\succ 0$ such that

Then for each outer iteration $s$, we propose $L$ number of inner iterations. This results in a nested block Jacobi method (NBJM) given in Algorithm 2.

Note that for the proposed splitting (17) of the positive definite block tri-diagonal matrix $\tilde{\Psi}$, condition $\rho(\tilde{\Phi}^{-1}\tilde{\Omega})<1$ always holds [20, Lem. 4.7 with Thm. 4.18].