Scalable computation for optimal control of
cascade systems with constraints

arXiv:1702.06375v2 [math.OC] 12 Oct 2017

Michael Cantoni, Farhad Farokhi, Eric Kerrigan, Iman Shames

∗

March 14, 2018

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.

1

Introduction

The application of model predictive control involves solving finite-horizon optimal control problems
in a receding horizon fashion [1–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–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–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
∗ Cantoni,

Farokhi and Shames are with the Department of Electrical and Electronic Engineering, The University
of Melbourne, Australia (Email: {cantoni,farhad.farokhi,iman.shames}@unimelb.edu.au). Kerrigan is with the
Departments of Electrical and Electronic Engineering and Aeronautics, Imperial College London, U.K. (Email:
e.kerrigan@imperial.ac.uk). This work is supported by the Australian Research Council (LP130100605) and a
McKenzie Fellowship.

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 ∈ N linear discrete-time dynamical sub-systems modelled by
xj (t + 1) = Aj (t)xj (t) + Bj (t)uj (t) + Ej (t)xj−1 (t),

(1)

given initial conditions xj (0) = ξj ∈ Rnj and model data Aj (t) ∈ Rnj ×nj , Bj (t) ∈ Rnj ×mj , and
Ej (t) ∈ Rnj ×nj−1 , with E1 (t) = 0 so that the spatial boundary value x0 (t) is effectively zero, for
j = 1, . . . , N and t = 0, . . . , T −1. The parameter T ∈ N is a specified time (or prediction) horizon.
The problem of interest is to set uj (t), for each sub-system index j = 1, . . . , N and sample time
(indexed by t = 0, . . . , T − 1), in order to minimize the separable penalty function
!


!
T
−1
N
X

 Qj (t) Sj (t)⊤ xj (t)
1X
⊤
⊤
⊤
xj (t)
uj (t)
+ xj (T ) Pj xj (T ) ,
(2)
J=
Sj (t) Rj (t) uj (t)
2
j=1

t=0

subject to separable inequality constraints
Mj (t)xj (t) + Lj (t)uj (t) ≤ cj (t) for j = 1, . . . , N, and t = 0, . . . , T,

(3)

where


Qj (t)
Sj (t)


Sj (t)⊤
 0,
Rj (t)

(4)

Mj (t) ∈ Rνj ×nj , Lj (t) ∈ Rνj ×mj and cj (t) ∈ Rνj , with 0 ≺ Rj (t) = Rj (t)⊤ ∈ Rmj ×mj , Mj (0) = 0,
Lj (T ) = 0, Sj (T ) = 0, and Qj (T ) = Pj = Pj⊤  0, are given for t = 0, . . . , T and j = 1, . . . , N .
Note that (4) implies Qj (t) = Qj (t)⊤  0 and Qj (t) − Sj (t)⊤ Rj (t)−1 Sj (t)  0, since Rj (t) ≻ 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 interiorpoint 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, . . . , N ,

ûj = uj (0)⊤

uj (T − 1)⊤

···

Then

⊤

∈ Rm j T

and x̂j


= xj (0)⊤

···

xj (T )⊤

⊤

∈ Rnj (T +1) .

−Âj x̂j + Êj x̂j−1 + B̂j ûj + Ĥj ξj = 0,

(5)

where


I


−Aj (0)


Âj =  0

 .
 ..
0

0
I
..
.
..

.
···

0
···

Bj (0) . . .


..
B̂j =  0
.

 .
.
..
 ..
0
···


···
..
.
..
.
..
.
0
···
..

.

..

.

0

···
..
..

.


0
.. 
.

..  ,
.


0



0


Ej (0)


Êj =  0

 .
 ..

.
0
−Aj (T − 1) I



0
I

..
0 

.
 

 .. 

..
. ,
,
and
Ĥ
=

j
.
 

.

 .. 

0
0
Bj (T − 1)

···
..
.
Ej (1)
..
.
···

···
..

···

.

..

.
0

..


0
.. 
.

..  ,
.

.. 
.

.
Ej (T − 1) 0

for j = 1, . . . , N . The boundary condition for (5) is effectively x̂0 = 0 as Ê1 = 0. Moreover, with
Q̂j = diag(Qj (0), . . . , Qj (T )) ∈ Rnj (T +1)×nj (T +1) ,
R̂j = diag(Rj (0), . . . , Rj (T − 1)) ∈ Rmj T ×mj T ,


Ŝj = diag(Sj (0), . . . , Sj (T )) 0 ∈ Rmj T ×nj (T +1) ,

M̂j = diag(Mj (0), . . . , Mj (T )) ∈ Rνj (T +1)×nj (T +1) ,

⊤
L̂j = diag(Lj (0), . . . , Lj (T − 1))⊤ 0 ∈ Rνj (T +1)×mj T ,


and ĉj = c⊤
j (0) · · ·
mulated as

⊤
c⊤
∈ Rνj (T +1) , the performance index and constraints can be reforj (T )
N

J=


1 X ⊤
⊤
x̂j Q̂j x̂j + û⊤
j R̂j ûj + 2ûj Ŝj x̂j
2 j=1

(6)

and
M̂j x̂j + L̂j ûj − ĉj ≤ 0 for j = 1, . . . , N,

(7)

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:
min

(û1 ,...,ûN )∈Rm1 T ×···×RmN T
(x̂1 ,...,x̂N )∈Rn1 (T +1) ×···×RnN (T +1)

(6)

subject to

(5) and (7).

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, . . . , N , where x̂0 = 0 (since Ê1 = 0), ÊN +1 = 0,
Λj = diag((λj )1 , . . . , (λj )νj ), Θj = diag((θj )1 , . . . , (θj )νj ), and 1 denotes a column vector of ones:
⊤
⊤
Q̂j x̂j + Ŝj⊤ ûj − Â⊤
j pj + M̂j λj + Êj+1 pj+1 = 0;

Ŝj x̂j + R̂j ûj + B̂j⊤ pj + L̂⊤
j λj = 0;
−Âj x̂j + Êj x̂j−1 + B̂j ûj + Ĥj ξj = 0;
M̂j x̂j + L̂j ûj − ĉj + θj = 0;


λ⊤
i

Λj Θj 1 = 0;
⊤
θj⊤ ≥ 0.


⊤

Given s[k] = s1 [k]⊤ · · · sN [k]⊤ , with sj [k] = x̂j [k]⊤ ûj [k]⊤
for j = 1, . . . , N , the iteration at Newton step k ∈ N is given by

pj [k]⊤

λj [k]⊤

θj [k]⊤

s[k + 1] = s[k] + α[k] · δ[k],

⊤

(8)

where α[k] > 0 is a sufficiently small step-size parameter selected online so that λj [k + 1] ≥ 0,

⊤
θj [k + 1] > 0, and δ[k] = δ1 [k]⊤ · · · δN [k]⊤ is the solution of the linearized KKT conditions

D1 [k] −Υ⊤
0
2

 −Υ2 D2 [k] −Υ⊤
3


..
 0
.
−Υ3

 .
..
..
 ..
.
.
0
···
0


 

δ1 [k]
ρ1 [k]
 δ [k]   ρ [k] 
 2   2 
 .   . 
 .   . 
.  =  . ,
0 

.   . 

 ..   .. 
DN −1 [k] −Υ⊤
N
ρN [k]
−ΥN
DN [k] δN [k]
...
..
.
..
.

0
..
.

(9)

and the following hold for j = 1, . . . , N , with ΥN +1 = 0 (since ÊN +1 = 0):


M̂j⊤
0
Q̂j Ŝj⊤ −Â⊤
j


L̂⊤
0 
B̂j⊤
R̂j
 Ŝj
j


Dj [k] = −Âj B̂j
0
0
0 ;


 M̂j L̂j
0
0
I 
0
0
0
Θj [k] Λj [k]


0
 0

Υj = 
−Êj
 0
0

0
0
0
0
0

0
0
0
0
0

0
0
0
0
0

h
ρj [k] = −Υj sj−1 [k] + Dj [k]sj [k] − Υ⊤
s
[k]
+
0⊤
j+1
j+1

(10)


0
0

0
;
0
0
0⊤

(11)

(Ĥj ξj )⊤

ĉ⊤
j

σj [k]⊤

i⊤

;

(12)

σj [k] = −Θj [k]λj [k] − Λj [k]θj [k] + Λj [k]Θj [k]1 − σ̄µ[k]1; Λj [k] = diag((λj [k])1 , . . . , (λj [k])νj );
PN
PN
Θj [k] = diag((θj [k])1 , . . . , (θj [k])νj ); σ̄ ∈ (0, 1) is a centring parameter; and µ[k] = ( j=1 (λj [k])⊤ θj [k])/( j=1 ν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 ∈ N, the unique
solution of (9) is given by the backward and forward recursions
(
(
ρN
for j = N
for j = 1
Σ−1
1 ρ̃1
ρ̃j =
,
and
δ
=
,
j
−1
−1
⊤
ρj + Υj+1 Σj+1 ρ̃j+1 for j = N − 1, . . . , 1
Σj (ρ̃j + Υj δj−1 ) for j = 2, . . . , N
respectively, where
Σj =

(

DN
−1
Dj − Υ ⊤
Υj+1
j+1 (Σj+1 )

for j = N
.
for j = N − 1, . . . , 1

(13)

In particular, each Σj ∈ Rpj ×pj , where pj = ((2nj + 2νj )(T + 1) + mj 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 Σj (c.f. Q̌j  0 in the proof of Lemma 1)
becomes a full matrix for j = N − 1, . . . , 1, despite the sparsity of Q̂j , R̂j , Ŝj , Âj , B̂j , M̂j , L̂j ,
Θj , and Λj . Therefore, inverting Σj in the recursions above incurs cost that scales cubically with
T , yielding an overall computation cost of order N T 3 .
The following result encapsulates a method for inverting the matrix Σj . This method involves
the inverses of block diagonal and block bi-diagonal matrices of order T × T , and the inverse of
one unstructured positive-definite (mj T ) × (mj T ) matrix; note that typically mj is smaller than
nj and νj . The computation cost is order T and T 3 , respectively, for an overall cost that scales
cubically with the time horizon. The matrix 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 Σj defined according to (13),
  
⊤
Q̂j − Êj+1
(Σ−1
X1
j+1 )33 Êj+1

X2  
Ŝj
  
=
X
Σj 
3
−
Âj
  
X4  
M̂j
X5
0

the unique solution of the linear equations
   
M̂j⊤ 0
Ŝj⊤ −Â⊤
X1
Y1
j

 Y2 
X
L̂⊤
0 
B̂j⊤
R̂j
2
j
 
 
 
X3 
B̂j
0
0
0 
 = Y3  ,



Y4 
L̂j
0
0
I  X4
X5
Y5
0
0
Θj Λj

(14)

given Y1 ∈ Rnj (T +1)×q , Y2 ∈ Rmj T ×q , Y3 ∈ Rνj (T +1)×q , Y4 ∈ Rνj (T +1)×q , Y5 ∈ Rνj (T +1)×q , can be
constructed as follows, for j = 1, . . . , N with ÊN +1 = 0:


   ⊤ −1
 
 
M̂j Θj Λj −M̂j⊤ Θ−1
Y1
X1
j
Y2  +  L̂⊤ Θ−1 Λj −L̂⊤ Θ−1  Y4 
X2  = Ψj Ωj Ψ⊤
j
j
j
j
j
Y5
Y3
X3
0
0
 
    
  
−Θ−1
Λj Θj−1
Y
X4
M̂j L̂j X1
j
+ 4
=
−
Y5
X2
X5
I
0
0
0

and

where


I

0
I
0

Ψj = −R̃j−1 S̃j
0

Z̃j
0
Ωj =  0
R̃j−1
−W̃j
0


0
,
−R̃j−1 B̃j⊤ Ã−⊤
j
Ã−⊤
j
⊤ 
−W̃j
0 ,
−W̃j Q̃j

(15)

(16)

Z̃j = Ã−1 B̃(R̃ + B̃ ⊤ Ã−⊤ Q̃Ã−1 B̃)−1 B̃ ⊤ Ã−⊤ ,

W̃j = I − Q̃j Z̃j ,
 "

⊤
Q̂j − Êj+1
(Σ−1
Q̄j S̃j⊤
j+1 )33 Êj+1
=
H̃j =
S̃j R̃j
Ŝj

(17)
(18)
Ŝj⊤
R̂j

#

"
#

M̂j⊤
−1
+
⊤ Θj Λj M̂j
L̂j

Q̃j = Q̄j − S̃j⊤ R̃j−1 S̃j  0, Ãj = Âj + B̂j R̃j−1 S̃j is non-singular and B̃j = B̂j .


L̂j  0,

(19)

Remark 2. The cost of computing Ã−1
j is order T because of the (lower) block bi-diagonal structure
−1
of Ãj . The (mj T ) × (mj T ) symmetric positive semi-definite matrix (R̃j + B̃j⊤ Ã−⊤
j Q̃j Ãj B̃j ) ≻ 0
is full, on the other hand. So the computation cost of inversion is order T 3 , whereby an approach
to computing Σj based on Lemma 2 scales as T 3 . An alternative approach based on inversion of
H̃j , assuming that it is non-singular by requiring
i
h
Q̂j Ŝj⊤
≻ 0,
Ŝj R̂j

for example, would also involve computation cost that scales as T 3 since the (nj (T +1))×(nj (T +1))
⊤
matrix Q̄j = Q̂j − Êj+1
(Σ−1
j+1 )33 Êj+1 in the 11-block of 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
Σ 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 ∈ {2, . . . , N − 1} needs to send Υj Σ−1
j ρ̃j to, and receive Υj δj−1
from, the agent for j − 1. Moreover, the processing agent for subsystem j ∈ {1, . . . , N − 1} only
needs access to 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

35
Direct solve (9) SPARSE OFF
Direct solve (9) SPARSE ON
Lemma 1 + Lemma 2
Method tailored to temporal structure

30

Computation time

25

20

15

10

5

0

0

5

10

15

20

25

30

35

40

Cascade length N (number of samples)

Figure 1: Computation time for 16 Newton iterations of the interior point method in the following
cases for varying N and fixed T = 5: (i) Solving (9) directly with sparsity exploitation disabled
(circle); (ii) Solving (9) directly with sparsity exploitation enabled (dashed); (iii) Solving (9)
via Lemma 1 and Lemma 2 (star); and (iv) Solving (9) with a method tailored to the temporal
structure (dot).
state-space model of the form (1) can be constructed for channels that operate under decentralized
distant-downstream control [11, 21]; 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 ∈ {1, . . . , N }, the following hold: n(j) = 4; m(j) = 1; and
ν(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

N = 10

45

Direct solve (9) SPARSE OFF
Direct solve (9) SPARSE ON
Lemma 1 + Lemma 2
Method tailored to temporal structure

40

35

Computation time

30

25

20

15

10

5

0

0

10

20

30

40

50

60

70

Time horizon T (samples)

Figure 2: Computation time for 16 Newton iterations of the interior point method in the following
cases for varying T and fixed N = 10: (i) Solving (9) directly with sparsity exploitation disabled
(circle); (ii) Solving (9) directly with sparsity exploitation enabled (dashed); (iii) Solving (9)
via Lemma 1 and Lemma 2 (star); and (iv) Solving (9) with a method tailored to the temporal
structure (dot).
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)\rho); (ii) Direct solution of (9)
with sparsity exploitation enabled (i.e., X=sparse(X); X\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 timehorizon 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.

A

A technical lemma and proofs

Lemma 3. Given Q = Q⊤ ∈ Rn×n , R = R⊤ ∈ Rm×m , S ∈ Rm×n , A ∈ Rn×n , B ∈ Rn×m ,
M ∈ Rν×n , L ∈ Rν×m , suppose that R ≻ 0, Q  0, Q − S ⊤ R−1 S  0, and Ã = A + BR−1 S is
non-singular. Let


Q S ⊤ −A⊤ M ⊤ 0
 S
R
B⊤
L⊤ 0 



D = −A B
0
0
0
,
M
L
0
0
I
0
0
0
Θ Λ
where Θ = diag(θ1 , . . . , θν ) ≻ 0 and Λ = diag(λ1 , . . . , λν )  0. Then D is non-singular. Moreover,
the 33-block of D−1 is given by
1

1

1

1

−Ã−⊤ Q̃ 2 (I + Q̃ 2 Ã−1 B̃ R̃−1 B̃ ⊤ Ã−⊤ Q̃ 2 )−1 Q̃ 2 Ã−1  0,
where


Q̄
S̃

 
S̃ ⊤
Q
=
S
R̃

  ⊤

S⊤
M
+
Θ−1 Λ M
R
L⊤


L  0,

Q̃ = Q̄ − S̃ ⊤ R̃−1 S̃  0 and B̃ = B.
Proof. First note that Θ ≻ 0 is invertible and that


0
Θ

I
Λ

As such, it follows that D is invertible

 
Q̄ S̃ ⊤ −A⊤
Q
 S̃
R̃
B⊤  =  S
−A
−A B
0

−1

=



−Θ−1 Λ
I


Θ−1
.
0

if and only if the Schur complement
  ⊤

−1 
S ⊤ −A⊤
M
0 
0 I
M
⊤ 
⊤


R
B
− L
0
Θ Λ
0
B
0
0
0

L 0
0 0



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 R̃ = R + L⊤ Θ−1 ΛL ≻ 0, since the diagonal matrix Θ−1 Λ  0 and R ≻ 0, and


Q̃ 0
−I
 0 R̃

(20)
0
−1
−1 ⊤ −⊤
−I 0 −Ã B̃ R̃ B̃ Ã




I
0
0
I
−S̃ ⊤ R̃−1
0
Q̄ S̃ ⊤ −A⊤
I
0   S̃
= 0
R̃
B ⊤  −R̃−1 S̃ I −R̃−1 B̃ ⊤ Ã−⊤  ,
0 −Ã−1 B̃ R̃−1 Ã−1
−A B
0
0
0
Ã−⊤
where Q̃ = Q̄ − S̃ ⊤ R̃−1 S̃  0 because
 
  ⊤

 
Q S⊤
M
Q̄ S̃ ⊤
M
−1
=
+
Θ
Λ
 0.
L
S R
L⊤
S̃ R̃
Therefore, D is non-singular if and only if
 

Q̃
−I
I
=
0
−I −Ã−1 B̃ R̃−1 B̃ ⊤ Ã−⊤

−Q̃
I



0
−I

−(I + Q̃Ã−1 B̃ R̃−1 B̃ ⊤ Ã−⊤ )
−Ã−1 B̃ R̃−1 B̃ ⊤ Ã−⊤



(21)

is non-singular, which is the case if and only if (I + Q̃Ã−1 B̃ R̃−1 B̃ ⊤ Ã−⊤ ) is non-singular, or
equivalently, if and only if −1 ∈
/ spec(Q̃Ã−1 B̃ R̃−1 B̃ ⊤ Ã−⊤ ). That later holds because
1

1

spec(Q̃Ã−1 B̃ R̃−1 B̃ ⊤ Ã−⊤ ) ∪ {0} = spec(Q̃ 2 Ã−1 B̃ R̃−1 B̃ ⊤ Ã−⊤ Q̃ 2 ) ∪ {0} ⊂ R≥0 ,
whereby D is invertible. In particular, the 22-block of the inverse of the left-hand side of (21)
1
1
1
1
is given by −(I + Q̃Ã−1 B̃ R̃−1 B̃ ⊤ Ã−⊤ )−1 Q̃ = −Q̃ 2 (I + Q̃ 2 Ã−1 B̃ R̃−1 B̃ ⊤ Ã−⊤ Q 2 )−1 Q̃ 2 , which is
congruent to the 33-block of D−1 via the transformation Ã−1 in view of (20), as claimed.

A.1

Proof of Lemma 1

Proof. Recall that R̂j ≻ 0, Q̂j  0, Q̂j − Ŝj⊤ R̂−1 Ŝj  0 and Ãj = Âj + B̂j R̂j−1 Ŝj is non-singular;
n.b., the structure of Ãj is the same as the block bi-diagonal structure of Âj with identity matrices
along the block diagonal. Now using Lemma 3, observe that Dj is invertible for j = 1, . . . , N .
Given the structure of Υj , the matrix Σj is the same as Dj except for the 11-block, which is Q̂j
−1
in the case of the latter and Q̌j = Q̂j − Êj⊤ (Σ−1
j+1 )33 Êj in the former. By Lemma 3, (DN )33  0,
so that (Σ−1
j+1 )33  0, whereby Q̌j  0, and thus, Σj is non-singular for j = N − 1 by Lemma 3
again. Continuing this argument inductively yields the invertibility of Σj for j = N − 2, . . . , 1. As
such, the specified recursions for Σj , ρ̃j and δj are all well defined.
−1
With ΣN = DN , solving the last block of (9) gives δN = Σ−1
N ρN + ΣN ΥN δN −1 . In turn,


D1


−Υ2


 0

 .
 ..
0

−Υ⊤
2

0

D2

−Υ⊤
3
..
.

−Υ3
..
.
···

..
0

.

...
..
.
..
.
DN −2
−ΥN −1

0
..
.



 

δ1
ρ1
 δ   ρ 
 2   2 
 .   . 
 .   . 
.  =  . ,
0 

.   . 

 ..   .. 
⊤

−ΥN −1
δN −1
ρ̃N −1
ΣN −1

where ΣN −1 and ρ̃N −1 are defined as in the statement of the lemma. Continuing in this fashion
yields the remaining expressions for Σj , ρ̃j and δj . The first block equation eventually becomes
Σ1 δ1k = ρ̃1 , and thus, δ1 = Σ−1
1 ρ̃1 .

A.2

Proof of Lemma 2

⊤
Proof. It is shown in the proof of Lemma 1 that −Êj+1
(Σ−1
j+1 )33 Êj+1  0. As such,



Q̂j
Ŝj

 
⊤
Ŝj⊤
−Êj+1
(Σ−1
j+1 )33 Êj+1
+
0
R̂j


0
 0,
0

and thus, Q̃j = Q̄j − S̃j⊤ R̃−1 S̃  0 as claimed. Also observe that Ãj is non-singular because it has
the same lower block bi-diagonal structure as Âj hwith identityimatrices along the block diagonal.

−1
−1
−1
Now recall that Θj ≻ 0, whereby Θ0j ΛIj
= −ΘI Λj Θ0 . Using this and the structure of

(14) yields


 
0
X4
=
Θj
X5

I
Λj

−1

and






Q̄j
X1
X2  =  S̃j
X3
−Âj



 
 Y4 − M̂j
Y5
0

L̂j
0

−1
−Â⊤
j
B̂j⊤ 


I

 0
0



S̃j⊤
R̃j
B̂j

0

 
 X1

0
0  
X2
=
Θj
0
X3

0
I
0

Q̄j
 S̃j
−Âj

S̃j⊤
R̃j
B̂j

−1 

#
  " ⊤
M̂j 0
0
0
−
⊤
0 
Θ
L̂j
0
j


I
0 0

In view of (20),


I
Λj



Y4
+ Θ−1
j Λj M̂j
Y5
 
Y1
−1 
Y2 
I

 

Λj
 
 Y3   .
 Y4 
Y5

−1 

I
0
0
−Â⊤
j
−1
−1
−⊤
B̂j⊤  = −R̃j S̃j I −R̃j B̃j⊤ Ãj 
−⊤
0
0
Ãj
0
−1 

Q̃j 0
−I
I


 0 R̃j
0
0
−1 ⊤ −⊤
−1
0
−I 0 −Ãj B̃j R̃j B̃j Ãj

−S̃j⊤ R̃j−1
I
−1
−Ãj B̃j R̃j−1

As such, the result follows by noting that


Q̃j
 0
−I

−1
0
−I

R̃j
0
−1 ⊤ −⊤
Ã
R̃
B̃
0 −Ã−1
B̃
j j
j
j
j
 −1
−1 ⊤ −⊤
0
Ãj B̃j R̃j B̃j Ãj W̃j
=
0
R̃−1
−W̃j

0


Ãj−1 B̃j R̃j−1 B̃j⊤ Ã−⊤
j W̃j Q̃j − I
,
0
−W̃j Q̃j

 
 X1
L̂j
X2

(22)


0
0 .
Ã−1
j

(23)

−1 ⊤ −⊤ −1
. Indeed, the matrix in (23) is precisely Ωj in (16).
where W̃j = (I + Q̃j Ã−1
j B̃j R̃j B̃j Ãj )
Application of the Sherman-Morrison-Woodbury matrix inversion lemma gives the equivalent ex−1 ⊤ −⊤
pression, including the expression (18) for W̃j , and (17) for Z̃j = Ã−1
j B̃j R̃j B̃j Ãj W̃j .

References
[1] C. Garcia, D. Prett, and M. Morari, “Model predictive control: Theory and practice – A
survey,” Automatica, vol. 25, no. 3, pp. 335–348, 1989.
[2] D. Q. Mayne, J. B. Rawlings, C. V. Rao, and P. O. M. Scokaert, “Constrained model predictive
control: Stability and optimality,” Automatica, vol. 36, no. 6, pp. 789–814, 2000.

[3] J. Maciejowski, Predictive Control: With Constraints. Prentice-Hall, 2002.
[4] J. B. Rawlings and D. Q. Mayne, Model Predictive Control: Theory and Design. Nob Hill
Pub., 2009.
[5] S. J. Wright, “Interior point methods for optimal control of discrete time systems,” Journal
of Optimization Theory and Applications, vol. 77, no. 1, pp. 161–187, 1993.
[6] C. V. Rao, S. J. Wright, and J. B. Rawlings, “Application of interior-point methods to model
predictive control,” Journal of Optimization Theory and Applications, vol. 99, no. 3, pp. 723–
757, 1998.
[7] M. Diehl, H. J. Ferreau, and N. Haverbeke, “Efficient numerical methods for nonlinear
MPC and moving horizon estimation,” in Nonlinear Model Predictive Control, pp. 391–417,
Springer, 2009.
[8] Y. Wang and S. Boyd, “Fast model predictive control using online optimization,” IEEE
Transactions on Control Systems Technology, vol. 18, no. 2, pp. 267–278, 2010.
[9] A. Domahidi, A. U. Zgraggen, M. N. Zeilinger, M. Morari, and C. N. Jones, “Efficient interior
point methods for multistage problems arising in receding horizon control,” in Proceedings
51st IEEE Conference on Decision and Control (CDC), pp. 668–674, IEEE, 2012.
[10] Y. Li, M. Cantoni, and E. Weyer, “On water-level error propagation in controlled irrigation
channels,” in Proceedings of the 44th IEEE Conference on Decision and Control, pp. 2101–
2106, IEEE, 2005.
[11] L. Soltanian and M. Cantoni, “Decentralized string-stability analysis for heterogeneous cascades subject to load-matching requirements,” Multidimensional Systems and Signal Processing, vol. 26, no. 4, pp. 985–999, 2015.
[12] V. Puig, G. Cembrano, J. Romera, J. Quevedo, B. Aznar, G. Ramón, and J. Cabot, “Predictive optimal control of sewer networks using CORAL tool: Application to Riera Blanca
catchment in Barcelona,” Water Science and Technology, vol. 60, no. 4, pp. 869–878, 2009.
[13] J. W. Labadie, “Optimal operation of multireservoir systems: State-of-the-art review,” Journal of Water Resources Planning and Management, vol. 130, no. 2, pp. 93–111, 2004.
[14] A. Şahin and M. Morari, “Decentralized model predictive control for a cascade of river power
plants,” in Intelligent Infrastructures, pp. 463–485, Springer, 2010.
[15] P. Seiler, A. Pant, and K. Hedrick, “Disturbance propagation in vehicle strings,” IEEE Transactions on Automatic Control, vol. 49, no. 10, pp. 1835–1841, 2004.
[16] S. M. Disney and D. R. Towill, “The effect of vendor managed inventory (VMI) dynamics on
the bullwhip effect in supply chains,” International Journal of Production Economics, vol. 85,
no. 2, pp. 199–215, 2003.
[17] S. Mehrotra, “On the implementation of a primal-dual interior point method,” SIAM Journal
on Optimization, vol. 2, no. 4, pp. 575–601, 1992.
[18] S. J. Wright, Primal-Dual Interior-Point Methods. SIAM, 1997.
[19] G. Meurant, “A review on the inverse of symmetric tridiagonal and block tridiagonal matrices,” SIAM Journal on Matrix Analysis and Applications, vol. 13, no. 3, pp. 707–728, 1992.
[20] R. Bevilacqua, B. Codenotti, and F. Romani, “Parallel solution of block tridiagonal linear
systems,” Linear Algebra and its Applications, vol. 104, pp. 39–57, 1988.

[21] M. Cantoni, E. Weyer, Y. Li, S. K. Ooi, I. Mareels, and M. Ryan, “Control of large-scale
irrigation networks,” Proceedings of the IEEE, vol. 95, no. 1, pp. 75–91, 2007.
[22] A. R. Neshastehriz, M. Cantoni, and I. Shames., “Water-level reference planning for automated irrigation channels via robust MPC,” in Proc. European Control Conference, pp. 1331–
1336, 2014.

