Scalable computation for optimal control of cascade systems with constraints

TL;DR

This paper introduces a scalable, distributed method for solving large-scale finite-horizon optimal control problems in cascade systems, with computation cost scaling linearly with the number of sub-systems.

Contribution

It develops a novel interior-point approach exploiting problem structure for efficient, distributed computation in cascade control systems.

Findings

Computation cost scales linearly with the number of sub-systems.

Method is suitable for distributed implementation with low communication overhead.

Example demonstrates application to automated irrigation channel dynamics.

Abstract

A method is devised for numerically solving a class of finite-horizon optimal control problems subject to cascade linear discrete-time dynamics. It is assumed that the linear state and input inequality constraints, and the quadratic measure of performance, are all separable with respect to the spatial dimension of the underlying cascade of sub-systems, as well as the temporal dimension of the dynamics. By virtue of this structure, the computation cost of an interior-point method for an equivalent quadratic programming formulation of the optimal control problem can be made to scale linearly with the number of sub-systems. However, the complexity of this approach grows cubically with the time horizon. As such, computational advantage becomes apparent in situations where the number of sub-systems is relatively large. In any case, the method is amenable to distributed computation with low communication overhead and only immediate upstream neighbour sharing of partial model data among processing agents. An example is presented to illustrate an application of the main results to model data for the cascade dynamics of an automated irrigation channel.

Equations

$$x_j(t+1)$$

$$J=\frac{1}{2}\sum _{j=1}^N\left(\left(\sum _{t=0}^{T-1}\left[\begin{array}{cc}{x}_j(t{)}^{\top }& u_j(t{)}^{\top }\end{array}\right]\left[\begin{array}{cc}{Q}_j(t)& S_j(t{)}^{\top }\\ S_j(t)& R_j(t)\end{array}\right]\left[\begin{array}{c}{x}_j(t)\\ u_j(t)\end{array}\right]\right)+x_j(T{)}^{\top }{P}_j{x}_j(T)\right),$$

$$M_j(t)x_j(t)$$

$$\left[\begin{array}{cc}{Q}_j(t)& S_j(t{)}^{\top }\\ S_j(t)& R_j(t)\end{array}\right]⪰0,$$

$${\hat{u}}_j$$

$$-{\hat{A}}_j{\hat{x}}_j+{\hat{E}}_j{\hat{x}}_{j-1}+{\hat{B}}_j{\hat{u}}_j+{\hat{H}}_j{\xi }_j=0,$$

$${\hat{A}}_j=\left[\begin{array}{ccccc}I& 0& \cdots & \cdots & 0\\ -A_j(0)& I& \ddots & & ⋮\\ 0& \ddots & \ddots & \ddots & ⋮\\ ⋮& \ddots & \ddots & \ddots & 0\\ 0& \cdots & 0& -A_j(T-1)& I\end{array}\right],{\hat{E}}_j=\left[\begin{array}{ccccc}0& \cdots & \cdots & \cdots & 0\\ E_j(0)& \ddots & & & ⋮\\ 0& E_j(1)& \ddots & & ⋮\\ ⋮& \ddots & \ddots & \ddots & ⋮\\ 0& \cdots & 0& E_j(T-1)& 0\end{array}\right],$$

$${\hat{B}}_j=\left[\begin{array}{cccc}0& \cdots & \cdots & 0\\ B_j(0)& \ddots & & ⋮\\ 0& \ddots & \ddots & ⋮\\ ⋮& \ddots & \ddots & 0\\ 0& \cdots & 0& B_j(T-1)\end{array}\right],\text{and}{\hat{H}}_j=\left[\begin{array}{c}I\\ 0\\ ⋮\\ ⋮\\ 0\end{array}\right],$$

$${\hat{Q}}_j$$

$${\hat{R}}_j$$

$${\hat{S}}_j$$

$${\hat{M}}_j$$

$${\hat{L}}_j$$

$$J=\frac{1}{2}\sum _{j=1}^N\left({\hat{x}}_j^{\top }{\hat{Q}}_j{\hat{x}}_j+{\hat{u}}_j^{\top }{\hat{R}}_j{\hat{u}}_j+2{\hat{u}}_j^{\top }{\hat{S}}_j{\hat{x}}_j\right)$$

$${\hat{M}}_j{\hat{x}}_j+{\hat{L}}_j{\hat{u}}_j-{\hat{c}}_j\le 0\text{ for }j=1,\dots ,N,$$

$$\underset{\begin{array}{c}({\hat{u}}_1,\dots ,{\hat{u}}_N)\in {\mathbb{R}}^{m_{1}T}\times \cdots \times {\mathbb{R}}^{m_{N}T}\\ ({\hat{x}}_1,\dots ,{\hat{x}}_N)\in {\mathbb{R}}^{n_1(T+1)}\times \cdots \times {\mathbb{R}}^{n_N(T+1)}\end{array}}{\min }(\text{ref}eq:cost2)\text{ subject to }(\text{ref}eq:spacemodel)\text{ and }(\text{ref}eq:constraints2).$$

$${\hat{Q}}_j{\hat{x}}_j+{\hat{S}}_j^{\top }{\hat{u}}_j-{\hat{A}}_j^{\top }{p}_j+{\hat{M}}_j^{\top }{\lambda }_j+{\hat{E}}_{j+1}^{\top }{p}_{j+1}$$

$${\hat{S}}_j{\hat{x}}_j+{\hat{R}}_j{\hat{u}}_j+{\hat{B}}_j^{\top }{p}_j+{\hat{L}}_j^{\top }{\lambda }_j$$

$$-{\hat{A}}_j{\hat{x}}_j+{\hat{E}}_j{\hat{x}}_{j-1}+{\hat{B}}_j{\hat{u}}_j+{\hat{H}}_j{\xi }_j$$

$${\hat{M}}_j{\hat{x}}_j+{\hat{L}}_j{\hat{u}}_j-{\hat{c}}_j+{\theta }_j$$

$${\Lambda }_j{\Theta }_{j}1$$

$${\left[\begin{array}{cc}{\lambda }_i^{\top }& {\theta }_j^{\top }\end{array}\right]}^{\top }$$

$$s[k+1]=s[k]+\alpha [k]\cdot \delta [k],$$

$$\left[\begin{array}{ccccc}{D}_1[k]& -{{\rm Y}}_2^{\top }& 0& \dots & 0\\ -{{\rm Y}}_2& D_2[k]& -{{\rm Y}}_3^{\top }& \ddots & ⋮\\ 0& -{{\rm Y}}_3& \ddots & \ddots & 0\\ ⋮& \ddots & \ddots & D_{N-1}[k]& -{{\rm Y}}_N^{\top }\\ 0& \cdots & 0& -{{\rm Y}}_N& D_N[k]\end{array}\right]\left[\begin{array}{c}{\delta }_1[k]\\ {\delta }_2[k]\\ ⋮\\ ⋮\\ {\delta }_N[k]\end{array}\right]=\left[\begin{array}{c}{\rho }_1[k]\\ {\rho }_2[k]\\ ⋮\\ ⋮\\ {\rho }_N[k]\end{array}\right],$$

$$D_j[k]$$

$${{\rm Y}}_j$$

$${\rho }_j[k]=-{{\rm Y}}_j{s}_{j-1}[k]+D_j[k]s_j[k]-{{\rm Y}}_{j+1}^{\top }{s}_{j+1}[k]+{\left[\begin{array}{ccccc}{0}^{\top }& 0^{\top }& ({\hat{H}}_j{\xi }_j{)}^{\top }& {\hat{c}}_j^{\top }& {\sigma }_j[k{]}^{\top }\end{array}\right]}^{\top };$$

$${\tilde{\rho }}_j=\{\begin{array}{ll}{\rho }_N& \text{for }j=N\\ {\rho }_j+{{\rm Y}}_{j+1}^{\top }{\Sigma }_{j+1}^{-1}{\tilde{\rho }}_{j+1}& \text{for }j=N-1,\dots ,1\end{array},\text{and }{\delta }_j=\{\begin{array}{ll}{\Sigma }_1^{-1}{\tilde{\rho }}_1& \text{for }j=1\\ {\Sigma }_j^{-1}({\tilde{\rho }}_j+{{\rm Y}}_j{\delta }_{j-1})& \text{for }j=2,\dots ,N\end{array},$$

$${\Sigma }_j$$

$${\Sigma }_j\left[\begin{array}{c}{X}_1\\ X_2\\ X_3\\ X_4\\ X_5\end{array}\right]=\left[\begin{array}{ccccc}{\hat{Q}}_j-{\hat{E}}_{j+1}^{\top }({\Sigma }_{j+1}^{-1}{)}_{33}{\hat{E}}_{j+1}& {\hat{S}}_j^{\top }& -{\hat{A}}_j^{\top }& {\hat{M}}_j^{\top }& 0\\ {\hat{S}}_j& {\hat{R}}_j& {\hat{B}}_j^{\top }& {\hat{L}}_j^{\top }& 0\\ -{\hat{A}}_j& {\hat{B}}_j& 0& 0& 0\\ {\hat{M}}_j& {\hat{L}}_j& 0& 0& I\\ 0& 0& 0& {\Theta }_j& {\Lambda }_j\end{array}\right]\left[\begin{array}{c}{X}_1\\ X_2\\ X_3\\ X_4\\ X_5\end{array}\right]=\left[\begin{array}{c}{Y}_1\\ Y_2\\ Y_3\\ Y_4\\ Y_5\end{array}\right],$$

$$\left[\begin{array}{c}{X}_1\\ X_2\\ X_3\end{array}\right]$$

$$\left[\begin{array}{c}{X}_4\\ X_5\end{array}\right]$$

$${\Psi }_j$$

$${\Omega }_j$$

$${\tilde{Z}}_j$$

$${\tilde{W}}_j$$

$${\tilde{H}}_j$$

$$\left[\begin{array}{cc}{\hat{Q}}_j& {\hat{S}}_j^{\top }\\ {\hat{S}}_j& {\hat{R}}_j\end{array}\right]\succ 0,$$

Full text

1 Introduction

The application of model predictive control involves solving finite-horizon optimal control problems in a receding horizon fashion [1, 2, 3, 4]. When the penalty function used to quantify performance and the inequality constraints on the system states and inputs all separate along the prediction horizon, additional structure in the equality constraints that encode the system dynamics can be exploited to devise efficient methods for computing the solutions. In particular, methods with computation costs that scale linearly with the time horizon and cubically with the number of states and inputs are well-known [5, 6, 7, 8, 9]. The cubic scaling of these methods in the spatial dimension of the problem, however, can be a limiting factor within the context of controlling large-scale interconnections of sub-systems. In this paper, interconnection structure is exploited over temporal structure to devise a more scalable method for problems with cascade dynamics in particular; i.e., when the system to control is the series interconnection of numerous sub-systems, each with a control input, an interconnection input, and an interconnection output. Such models arise in the study of irrigation and drainage networks [10, 11, 12], mutli-reservoir and hydro-power systems [13, 14], vehicle platoons [15], and supply chain management [16], for example.

The proposed method for solving finite-horizon optimal control problems with cascade dynamics is closely related to the interior-point method developed in [5, 6]. Interchanging the roles of the temporal and spatial dimensions of such problems yields linear scaling of computation cost with the number of sub-systems along the cascade. However, the complexity grows cubically with the time horizon, despite the causal flow of information in the temporal dimension. The development illuminates the difficulty of overcoming such cubic dependence. Computational advantage over methods that exploit temporal structure, rather than the spatial structure exploited here, arises when the length of the cascade is relatively large compared to the prediction horizon. In any case, the method is amenable to distributed computation over a linear network of processing agents, one for each sub-system, with limited neighbour-to-neighbour communication, and only partial sharing of model information between neighbouring agents.

The paper is organized as follows. The class of finite-horizon optimal control problems studied is defined in Section 2.1. The formulation of an equivalent quadratic program (QP) is given in Section 2.2. A scalable interior point method for computing an optimal solution of the structured QP is developed in Section 3. Proofs are deferred to the Appendix. Finally, a numerical example based on model data for an automated irrigation channel is presented Section 4. Some concluding remarks are provided in Section 5.

2 Problem Formulation

A class of finite-horizon optimal control problems with cascade dynamics is defined in this section. In addition to the directed interconnection structure of the underlying cascade of sub-systems, a defining characteristic of this class is the separability of the state and input inequality constraints and performance index across both the spatial and temporal dimensions of the system dynamics. An equivalent QP with computationally favorable structure is formulated in Section 2.2.

2.1 Constrained finite-horizon optimal control

Consider the cascade of $N\in\mathbb{N}$ linear discrete-time dynamical sub-systems modelled by

[TABLE]

given initial conditions $x_{j}(0)=\xi_{j}\in\mathbb{R}^{n_{j}}$ and model data $A_{j}(t)\in\mathbb{R}^{n_{j}\times n_{j}}$, $B_{j}(t)\in\mathbb{R}^{n_{j}\times m_{j}}$, and $E_{j}(t)\in\mathbb{R}^{n_{j}\times n_{j-1}}$, with $E_{1}(t)=0$ so that the spatial boundary value $x_{0}(t)$ is effectively zero, for $j=1,\ldots,N$ and $t=0,\ldots,T-1$. The parameter $T\in\mathbb{N}$ is a specified time (or prediction) horizon. The problem of interest is to set $u_{j}(t)$, for each sub-system index $j=1,\ldots,N$ and sample time (indexed by $t=0,\ldots,T-1$), in order to minimize the separable penalty function

[TABLE]

subject to separable inequality constraints

[TABLE]

where

[TABLE]

$M_{j}(t)\in\mathbb{R}^{\nu_{j}\times n_{j}}$, $L_{j}(t)\in\mathbb{R}^{\nu_{j}\times m_{j}}$ and $c_{j}(t)\in\mathbb{R}^{\nu_{j}}$, with $0\prec R_{j}(t)=R_{j}(t)^{\top}\in\mathbb{R}^{m_{j}\times m_{j}}$, $M_{j}(0)=0$, $L_{j}(T)=0$, $S_{j}(T)=0$, and $Q_{j}(T)=P_{j}=P_{j}^{\top}\succeq 0$, are given for $t=0,\ldots,T$ and $j=1,\ldots,N$. Note that (4) implies $Q_{j}(t)=Q_{j}(t)^{\top}\succeq 0$ and $Q_{j}(t)-S_{j}(t)^{\top}R_{j}(t)^{-1}S_{j}(t)\succeq 0$, since $R_{j}(t)\succ 0$.

It is possible to reformulate the optimal control problem defined above in a number of ways that result in standard QPs. Following the style of QP reformulation in [5] leads to an interior-point method involving the solution of linear algebra problems with favourable block tridiagonal structure. This is exploited in Section 3 to devise an algorithm with per-iteration computation cost that scales linearly with cascade length $N$.

2.2 A QP formulation

First note that the equality constraint corresponding to the dynamics (1) can be reformulated as follows. Define, for $j=1,\ldots,N$,

[TABLE]

Then

[TABLE]

where

[TABLE]

for $j=1,\ldots,N$. The boundary condition for (5) is effectively $\hat{x}_{0}=0$ as $\hat{E}_{1}=0$. Moreover, with

[TABLE]

and $\hat{c}_{j}=\begin{bmatrix}c_{j}^{\top}(0)&\cdots&c_{j}^{\top}(T)\end{bmatrix}^{\top}\in\mathbb{R}^{\nu_{j}(T+1)}$, the performance index and constraints can be reformulated as

[TABLE]

and

[TABLE]

respectively. To summarize, the problem of minimizing (2) subject to the equality constraints (1), and inequality constraints (3), is now in the form of the following standard QP:

[TABLE]

3 Developing an interior-point method that scales linearly

with cascade length

A primal-dual interior-point method for solving a QP involves the application of Newton’s iterative method to solve the Karush–Kuhn–Tucker (KKT) conditions [17, 18]. For the case at hand, the KKT conditions take the following form, for $j=1,\ldots,N$, where $\hat{x}_{0}=0$ (since $\hat{E}_{1}=0$), $\hat{E}_{N+1}=0$, $\Lambda_{j}=\mathrm{diag}((\lambda_{j})_{1},\ldots,(\lambda_{j})_{\nu_{j}})$, $\Theta_{j}=\mathrm{diag}((\theta_{j})_{1},\ldots,(\theta_{j})_{\nu_{j}})$, and $\mathbf{1}$ denotes a column vector of ones:

[TABLE]

Given $s[k]=\begin{bmatrix}s_{1}[k]^{\top}&\cdots&s_{N}[k]^{\top}\end{bmatrix}^{\top}$, with $s_{j}[k]=\begin{bmatrix}\hat{x}_{j}[k]^{\top}&\hat{u}_{j}[k]^{\top}&p_{j}[k]^{\top}&\lambda_{j}[k]^{\top}&\theta_{j}[k]^{\top}\end{bmatrix}^{\top}$ for $j=1,\ldots,N$, the iteration at Newton step $k\in\mathbb{N}$ is given by

[TABLE]

where $\alpha[k]>0$ is a sufficiently small step-size parameter selected online so that $\lambda_{j}[k+1]\geq 0$, $\theta_{j}[k+1]>0$, and $\delta[k]=\begin{bmatrix}\delta_{1}[k]^{\top}&\cdots&\delta_{N}[k]^{\top}\end{bmatrix}^{\top}$ is the solution of the linearized KKT conditions

[TABLE]

and the following hold for $j=1,\ldots,N$, with $\Upsilon_{N+1}=0$ (since $\hat{E}_{N+1}=0$):

[TABLE]

[TABLE]

[TABLE]

$\sigma_{j}[k]=-\Theta_{j}[k]\lambda_{j}[k]-\Lambda_{j}[k]\theta_{j}[k]+\Lambda_{j}[k]\Theta_{j}[k]\mathbf{1}-\bar{\sigma}\mu[k]\mathbf{1}$; $\Lambda_{j}[k]=\mathrm{diag}((\lambda_{j}[k])_{1},\ldots,(\lambda_{j}[k])_{\nu_{j}})$; $\Theta_{j}[k]=\mathrm{diag}((\theta_{j}[k])_{1},\ldots,(\theta_{j}[k])_{\nu_{j}})$; $\bar{\sigma}\in(0,1)$ is a centring parameter; and $\mu[k]=(\sum_{j=1}^{N}(\lambda_{j}[k])^{\top}\theta_{j}[k])/(\sum_{j=1}^{N}\nu_{j})$ is the duality gap. Given an appropriate initialization $s[1]$, the linear equation (9) has a unique solution for each iteration (8), as seen below. Construction of this solution is facilitated by the block tridiagonal structure; e.g., see [19] and [20]. Indeed, the computation cost of solving the set of equations (9) at each iteration can be made to scale linearly with $N$; ignoring structure would incur order $N^{3}$ complexity. Proofs of the following results, which underpin this assertion, can be found in the Appendix.

Lemma 1**.**

Dropping explicit dependence on the algorithm iteration index $k\in\mathbb{N}$, the unique solution of (9) is given by the backward and forward recursions

[TABLE]

respectively, where

[TABLE]

In particular, each $\Sigma_{j}\in\mathbb{R}^{p_{j}\times p_{j}}$, where $p_{j}=((2n_{j}+2\nu_{j})(T+1)+m_{j}T)$, is non-singular as required.

Remark 1**.**

Using Lemma 1 to solve (9) incurs a computation cost that scales linearly with the number of sub-systems $N$. Unfortunately, the $11$-block of $\Sigma_{j}$ (c.f. $\check{Q}_{j}\succeq 0$ in the proof of Lemma 1) becomes a full matrix for $j=N-1,\ldots,1$, despite the sparsity of $\hat{Q}_{j}$, $\hat{R}_{j}$, $\hat{S}_{j}$, $\hat{A}_{j}$, $\hat{B}_{j}$, $\hat{M}_{j}$, $\hat{L}_{j}$, $\Theta_{j}$, and $\Lambda_{j}$. Therefore, inverting $\Sigma_{j}$ in the recursions above incurs cost that scales cubically with $T$, yielding an overall computation cost of order $NT^{3}$.

The following result encapsulates a method for inverting the matrix $\Sigma_{j}$. This method involves the inverses of block diagonal and block bi-diagonal matrices of order $T\times T$, and the inverse of one unstructured positive-definite $(m_{j}T)\times(m_{j}T)$ matrix; note that typically $m_{j}$ is smaller than $n_{j}$ and $\nu_{j}$. The computation cost is order $T$ and $T^{3}$, respectively, for an overall cost that scales cubically with the time horizon. The matrix $\hat{Q}_{j}$ is not required to be non-singular, as needed to follow steps used to invert similarly structured matrices in [8, 9], for example. Furthermore, the method is amenable to distributed implementation with low communication overhead, as also discussed in more detail subsequently.

Lemma 2**.**

With $\Sigma_{j}$ defined according to (13), the unique solution of the linear equations

[TABLE]

given $Y_{1}\in\mathbb{R}^{n_{j}(T+1)\times q}$, $Y_{2}\in\mathbb{R}^{m_{j}T\times q}$, $Y_{3}\in\mathbb{R}^{\nu_{j}(T+1)\times q}$, $Y_{4}\in\mathbb{R}^{\nu_{j}(T+1)\times q}$, $Y_{5}\in\mathbb{R}^{\nu_{j}(T+1)\times q}$, can be constructed as follows, for $j=1,\ldots,N$ with $\hat{E}_{N+1}=0$:

[TABLE]

*where *

[TABLE]

$\tilde{Q}_{j}=\bar{Q}_{j}-\tilde{S}_{j}^{\top}\tilde{R}_{j}^{-1}\tilde{S}_{j}\succeq 0$, $\tilde{A}_{j}=\hat{A}_{j}+\hat{B}_{j}\tilde{R}_{j}^{-1}\tilde{S}_{j}$ is non-singular and $\tilde{B}_{j}=\hat{B}_{j}$.

Remark 2**.**

The cost of computing $\tilde{A}_{j}^{-1}$ is order $T$ because of the (lower) block bi-diagonal structure of $\tilde{A}_{j}$. The $(m_{j}T)\times(m_{j}T)$ symmetric positive semi-definite matrix $(\tilde{R}_{j}+\tilde{B}_{j}^{\top}\tilde{A}_{j}^{-\top}\tilde{Q}_{j}\tilde{A}_{j}^{-1}\tilde{B}_{j})\succ 0$ is full, on the other hand. So the computation cost of inversion is order $T^{3}$, whereby an approach to computing $\Sigma_{j}$ based on Lemma 2 scales as $T^{3}$. An alternative approach based on inversion of $\tilde{H}_{j}$, assuming that it is non-singular by requiring

[TABLE]

*for example, would also involve computation cost that scales as $T^{3}$ since the $(n_{j}(T+1))\times(n_{j}(T+1))$ matrix $\bar{Q}_{j}=\hat{Q}_{j}-\hat{E}_{j+1}^{\top}(\Sigma_{j+1}^{-1})_{33}\hat{E}_{j+1}$ in the $11$-block of $\tilde{H}_{j}$ is full for $j<N$. By constrast, as seen in [8, 9], for example, taking such an approach can be fruitful within the context of exploiting temporal structure in optimal control problems, since the $11$-block (i.e., $Q$ block) of the corresponding $\Sigma$ matrix is sparse in this case. *

Finally, it is of note that the calculations required to implement the solution of (9) according to Lemma 1 can be distributed among processing agents, one for each sub-system, connected in a linear network that mirrors the underlying cascade of sub-system models. Indeed, the agent associated with sub-system $j\in\{2,\ldots,N-1\}$ needs to send $\Upsilon_{j}\Sigma_{j}^{-1}\tilde{\rho}_{j}$ to, and receive $\Upsilon_{j}\delta_{j-1}$ from, the agent for $j-1$. Moreover, the processing agent for subsystem $j\in\{1,\ldots,N-1\}$ only needs access to $\hat{E}_{j+1}$ (i.e., the influence of $j$ on $j+1$) and no further model, performance index or constraint data from other sub-systems. The processing agent for $j=1$ need only communicate with the agent for $j=2$ and the processing agent for $j=N$ need only communicate with the agent for $j=N-1$. For each Newton iteration of the interior-point method, information is communicated in a sequence of steps, once up the cascade, and then once down the cascade, without further iteration. Calculation of the the duality gap and step size updates can be determined by a similar pattern of up and then down exchange of information. Agent pipelining by the provision of buffering in the communication could be exploited to improve throughput, without of course improving latency, which may degrade slightly. Of course, it is necessary to process a sufficient number of Newton updates in the interior-point method. One of the appealing features of interior point methods is the typically small number of iterations (e.g., ten to fifteen) needed to reach a good solution to the QP. As such, the overall communication overhead of a distributed version of the method is low. The localization of model and other problem data can also be considered advantageous from security and privacy perspectives.

4 Example

Numerical results are obtained by applying the preceding developments to model data for an automated irrigation channel, within the context of a water-level reference planning problem. A state-space model of the form (1) can be constructed for channels that operate under decentralized distant-downstream control [21, 11]; each sub-system corresponds to a stretch of channel between adjacent flow regulators called a pool. An optimal control problem can be formulated to determine the water-level reference input for each pool, across a planning horizon for which a load forecast may be known, subject to hard constraints on the water-levels and flows. For this example, the model data of the pools, the distributed controllers, the discretization sample-period, and constraint levels are all borrowed from [22]; the many pool channels considered here are constructed by concatenating sections of the channel considered there. In the model for each sub-system (i.e., pool) there is one control input, the water-level reference, and four states, including two for the water-level dynamics and two for the PI-type feedback controller that sets the upstream inflow on the basis of the measured downstream water-level error. Box type constraints on the water-level and controller flow output states are to be satisfied in addition to the water-level reference input constraints. For each sub-system $j\in\{1,\ldots,N\}$, the following hold: $n(j)=4$; $m(j)=1$; and $\nu(j)=6$. The other model parameters (e.g., entries of state-space matrices) are not uniform along the channel.

Figure 1 shows the MATLAB 2014b computation time (in seconds with forced single computation thread) of exactly sixteen iterations of the interior point method described above for different ways of solving (9). The number of sub-systems $N$ is varied from $1$ to $40$, with fixed time-horizon $T=5$. In all cases the duality gap is less than $10^{-3}$ after sixteen iterations. The following cases are considered: $(i)$ Direct solution of (9) with sparsity exploitation disabled by perturbing the block tridiagonal matrix $X$ on the left-hand side (i.e., (X+eps)$\backslash$rho); $(ii)$ Direct solution of (9) with sparsity exploitation enabled (i.e., X=sparse(X); X$\backslash$rho); $(iii)$ Solution of (9) via Lemma 1 and Lemma 2; and $(iv)$ Solution of (9) via a method tailored to exploit the temporal structure also present in $X$, along the lines of Lemma 1. As expected, the approach that does not exploit structure incurs a computation time that grows as $N^{3}$. The use of Lemma 1 with Lemma 2, on the other hand, scales linearly with $N$. Moreover, the performance achieved is as good as enabling MATLAB to exploit structure when solving (9) directly, which of course requires the MATLAB environment. Also note that a method tailored to the temporal structure of $X$ does not scale linearly. By contrast, Figure 2 shows the computation time for fixed $N=10$ and varying time-horizon $T$. As expected, the approaches in cases $(i)$ and $(iii)$ scale as $T^{3}$. This is similar to the cubic scaling of computation cost with the dimension of the state in methods that are tailored to exploit the temporal structure of optimal control problems; e.g., [8, 9]. Direct solution of (9) with sparsity exploitation enabled appears to asymptotically scale linearly with $T$, in a way that is consistent with the method that is tailored to exploit temporal structure.

5 Conclusion

The main contribution is a scalable interior-point method for computing the solution of constrained discrete-time optimal control problems with cascade dynamics, over a finite prediction horizon. By exploiting the spatial structure arising from the cascade dynamics, the computation cost of each step scales linearly with the number of sub-systems along the cascade. By constrast, the method exhibits cubic growth of computation time as the prediction horizon increases. Direct application of standard methods, typically tailored to exploit the temporal structure of optimal control problems in order to achieve linear scaling with the time horizon, would yield complexity that grows as the cube of the number of system states, and thus, the number of sub-systems. The main developments are illustrated by numerical example on model data for an automated irrigation channel. A topic of ongoing research pertains to extending the main ideas to exploit directed and undirected spatial propagation of information in tree networks of dynamical systems. Another concerns the design of custom hardware for distributed algorithm implementation.

Appendix A A technical lemma and proofs

Lemma 3**.**

Given $Q=Q^{\top}\in\mathbb{R}^{n\times n}$, $R=R^{\top}\in\mathbb{R}^{m\times m}$, $S\in\mathbb{R}^{m\times n}$, $A\in\mathbb{R}^{n\times n}$, $B\in\mathbb{R}^{n\times m}$, $M\in\mathbb{R}^{\nu\times n}$, $L\in\mathbb{R}^{\nu\times m}$, suppose that $R\succ 0$, $Q\succeq 0$, $Q-S^{\top}R^{-1}S\succeq 0$, and $\tilde{A}=A+BR^{-1}S$ is non-singular. Let

[TABLE]

where $\Theta=\mathrm{diag}(\theta_{1},\ldots,\theta_{\nu})\succ 0$ and $\Lambda=\mathrm{diag}(\lambda_{1},\ldots,\lambda_{\nu})\succeq 0$. Then $D$ is non-singular. Moreover, the $33$-block of $D^{-1}$ is given by

[TABLE]

where

[TABLE]

$\tilde{Q}=\bar{Q}-\tilde{S}^{\top}\tilde{R}^{-1}\tilde{S}\succeq 0$* and $\tilde{B}=B$.*

Proof.

First note that $\Theta\succ 0$ is invertible and that

[TABLE]

As such, it follows that $D$ is invertible if and only if the Schur complement

[TABLE]

is non-singular. Moreover, the inverse of this matrix, when it exists, is precisely the matrix comprising the first three block rows and columns of $D^{-1}$; as such, the $33$-blocks coincide.

Now note that $\tilde{R}=R+L^{\top}\Theta^{-1}\Lambda L\succ 0$, since the diagonal matrix $\Theta^{-1}\Lambda\succeq 0$ and $R\succ 0$, and

[TABLE]

where $\tilde{Q}=\bar{Q}-\tilde{S}^{\top}\tilde{R}^{-1}\tilde{S}\succeq 0$ because

[TABLE]

Therefore, $D$ is non-singular if and only if

[TABLE]

is non-singular, which is the case if and only if $(I+\tilde{Q}\tilde{A}^{-1}\tilde{B}\tilde{R}^{-1}\tilde{B}^{\top}\tilde{A}^{-\top})$ is non-singular, or equivalently, if and only if $-1\notin\mathrm{spec}(\tilde{Q}\tilde{A}^{-1}\tilde{B}\tilde{R}^{-1}\tilde{B}^{\top}\tilde{A}^{-\top})$. That later holds because

[TABLE]

whereby $D$ is invertible. In particular, the $22$-block of the inverse of the left-hand side of (21) is given by $-(I+\tilde{Q}\tilde{A}^{-1}\tilde{B}\tilde{R}^{-1}\tilde{B}^{\top}\tilde{A}^{-\top})^{-1}\tilde{Q}=-\tilde{Q}^{\frac{1}{2}}(I+\tilde{Q}^{\frac{1}{2}}\tilde{A}^{-1}\tilde{B}\tilde{R}^{-1}\tilde{B}^{\top}\tilde{A}^{-\top}Q^{\frac{1}{2}})^{-1}\tilde{Q}^{\frac{1}{2}}$, which is congruent to the $33$-block of $D^{-1}$ via the transformation $\tilde{A}^{-1}$ in view of (20), as claimed. ∎

A.1 Proof of Lemma 1

Proof.

Recall that $\hat{R}_{j}\succ 0$, $\hat{Q}_{j}\succeq 0$, $\hat{Q}_{j}-\hat{S}_{j}^{\top}\hat{R}^{-1}\hat{S}_{j}\succeq 0$ and $\tilde{A}_{j}=\hat{A}_{j}+\hat{B}_{j}\hat{R}_{j}^{-1}\hat{S}_{j}$ is non-singular; n.b., the structure of $\tilde{A}_{j}$ is the same as the block bi-diagonal structure of $\hat{A}_{j}$ with identity matrices along the block diagonal. Now using Lemma 3, observe that $D_{j}$ is invertible for $j=1,\ldots,N$. Given the structure of $\Upsilon_{j}$, the matrix $\Sigma_{j}$ is the same as $D_{j}$ except for the $11$-block, which is $\hat{Q}_{j}$ in the case of the latter and $\check{Q}_{j}=\hat{Q}_{j}-\hat{E}_{j}^{\top}(\Sigma_{j+1}^{-1})_{33}\hat{E}_{j}$ in the former. By Lemma 3, $(D_{N}^{-1})_{33}\preceq 0$, so that $(\Sigma_{j+1}^{-1})_{33}\preceq 0$, whereby $\check{Q}_{j}\succeq 0$, and thus, $\Sigma_{j}$ is non-singular for $j=N-1$ by Lemma 3 again. Continuing this argument inductively yields the invertibility of $\Sigma_{j}$ for $j=N-2,\ldots,1$. As such, the specified recursions for $\Sigma_{j}$, $\tilde{\rho}_{j}$ and $\delta_{j}$ are all well defined.

With $\Sigma_{N}=D_{N}$, solving the last block of (9) gives $\delta_{N}=\Sigma_{N}^{-1}\rho_{N}+\Sigma_{N}^{-1}\Upsilon_{N}\delta_{N-1}$. In turn,

[TABLE]

where $\Sigma_{N-1}$ and $\tilde{\rho}_{N-1}$ are defined as in the statement of the lemma. Continuing in this fashion yields the remaining expressions for $\Sigma_{j}$, $\tilde{\rho}_{j}$ and $\delta_{j}$. The first block equation eventually becomes $\Sigma_{1}\delta_{1}^{k}=\tilde{\rho}_{1}$, and thus, $\delta_{1}=\Sigma_{1}^{-1}\tilde{\rho}_{1}$. ∎

A.2 Proof of Lemma 2

Proof.

It is shown in the proof of Lemma 1 that $-\hat{E}_{j+1}^{\top}(\Sigma_{j+1}^{-1})_{33}\hat{E}_{j+1}\succeq 0$. As such,

[TABLE]

and thus, $\tilde{Q}_{j}=\bar{Q}_{j}-\tilde{S}_{j}^{\top}\tilde{R}^{-1}\tilde{S}\succeq 0$ as claimed. Also observe that $\tilde{A}_{j}$ is non-singular because it has the same lower block bi-diagonal structure as $\hat{A}_{j}$ with identity matrices along the block diagonal. Now recall that $\Theta_{j}\succ 0$, whereby $\left[\begin{smallmatrix}0&I\\ \Theta_{j}&\Lambda_{j}\end{smallmatrix}\right]^{-1}=\left[\begin{smallmatrix}-\Theta^{-1}\!\Lambda_{j}&~{}\Theta^{-1}\\ I&0\end{smallmatrix}\right]$. Using this and the structure of (14) yields

[TABLE]

and

[TABLE]

In view of (20),

[TABLE]

As such, the result follows by noting that

[TABLE]

where $\tilde{W}_{j}=(I+\tilde{Q}_{j}\tilde{A}_{j}^{-1}\tilde{B}_{j}\tilde{R}_{j}^{-1}\tilde{B}_{j}^{\top}\tilde{A}_{j}^{-\top})^{-1}$. Indeed, the matrix in (23) is precisely $\Omega_{j}$ in (16). Application of the Sherman-Morrison-Woodbury matrix inversion lemma gives the equivalent expression, including the expression (18) for $\tilde{W}_{j}$, and (17) for $\tilde{Z}_{j}=\tilde{A}_{j}^{-1}\tilde{B}_{j}\tilde{R}_{j}^{-1}\tilde{B}_{j}^{\top}\tilde{A}_{j}^{-\top}\tilde{W}_{j}$. ∎