A Complexity Analysis of Event-Triggered Model Predictive Control on Industrial Hardware

TL;DR

This paper evaluates the practical implementation of event-triggered model predictive control on industrial hardware, analyzing computational and bandwidth trade-offs through theoretical and empirical methods.

Contribution

It provides a detailed complexity analysis of different implementation variants, combining theoretical insights with real-world hardware experiments.

Findings

Bandwidth is minimized when transmitting only active constraints.

Computational cost scales nearly linearly with problem size in practice.

Implementation results confirm theoretical complexity analysis.

Abstract

We implement a recently proposed event-triggered networked MPC approach on industrial hardware to analyze its practical relevance. There exist several alternatives for such an implementation that differ with respect to the distribution of computational load between local and central nodes, and with respect to network bandwidth requirements. These alternatives have been analyzed theoretically before, but when implemented it becomes evident that their usefulness cannot be predicted based on theoretical considerations alone. It is the purpose of the present paper to account for both practical and theoretical aspects in determining which alternative is most appropriate for an implementation on industrial hardware. The smallest possible bandwidth is known to result for a variant in which only the active set of constraints is transmitted from the central to the local nodes. Since local nodes must determine the control law from the active set in this case, which requires matrix inversions, an unattractive computational cost results at first sight. Somewhat surprisingly, the computational cost scales practically linearly in the problem size when implemented. We confirm this result with a more detailed theoretical complexity analysis than given in previous papers. All results are illustrated with data obtained with an implementation on industrial hardware components.

Tables

Table 1. Table 1: Computational effort for some matrix operations

op.	$M+\tilde{M}$	$cM$	$MV$	$M^{-1}$
flops	$mn$	$mn$	$ml\left(2n-1\right)$	$\frac{2n^3+18n^2+10n}{3}$

Equations

$$\mathcal{P}=\{x\in {\mathbb{R}}^n\mid Tx\le d\}.$$

$$x(0)$$

$$x(N)$$

$$\underset{U}{\min }\frac{1}{2}{U}^{\prime }HU+x^{\prime }FU\text{s.t. }GU-Ex\le w$$

$$U^{\star}(x)=\left\{\begin{array}[]{cl}K_{1}x+b_{1}\text{, if }x\in\mathcal{P}_{1},\\ \vdots\\ K_{p}x+b_{p}\text{, if }x\in\mathcal{P}_{p},\end{array}\right.$$

$$\begin{array}{rl}\mathcal{A}& =\left\{i\in \left\{1,...,q\right\}|G_i{U}^{\star }(x)-E_{i}x=w_i\right\},\\ \mathcal{I}& =\left\{i\in \left\{1,...,q\right\}|G_i{U}^{\star }(x)-E_{i}x<w_i\right\}.\end{array}$$

$$\displaystyle\scalebox{0.93}{$K_{\mathcal{M}}^{\star}$}$$

$$\displaystyle\scalebox{0.93}{$b_{\mathcal{M}}^{\star}$}$$

$$\displaystyle\scalebox{0.93}{$T^{\star}$}$$

$$\displaystyle\scalebox{0.93}{$d^{\star}$}$$

$$u^{\star }(x)$$

$${\mathcal{P}}^{\star }$$

$$e(x)=\left\{\begin{array}[]{ll}0&\text{ if }x\in\mathcal{P}^{\star},\\ 1&\text{ otherwise.}\end{array}\right.$$

$$\mathcal{A}(x),$$

$$U^{\star }(x),$$

$$K_{\mathcal{M}}^{\star },b_{\mathcal{M}}^{\star },T^{\star },d^{\star }.$$

$$\gamma_{i}=\left\{\begin{array}[]{cl}0&\mbox{if }G_{i}U^{\star}+E_{i}x<w_{i}\\ 1&\mbox{if }G_{i}U^{\star}+E_{i}x=w_{i}\end{array}\right.,$$

$$\Phi =(G_{\mathcal{A}}{H}^{-1}(G_{\mathcal{A}}{)}^{\prime }{)}^{-1}$$

$$\frac{\lambda -2}{3}>\frac{n}{m}$$

$$\frac{\lambda -2}{3}>\frac{n}{m}+\frac{\lambda (q_{\mathcal{A}}(q_{\mathcal{A}}+1))}{6mN}.$$

$$\lambda mN>q.$$

$$\lambda >2+2\frac{n}{m}(1+\frac{1}{N}).$$

$$\lambda >2+3\frac{n}{m}$$

$$\frac{{\eta }_{\mathrm{inv}}(q_{\mathcal{A}})}{{\eta }_{\mathrm{mat}}(q_{\mathcal{A}})}\le \frac{18}{79}.$$

$${\eta }_{\mathrm{inv}}(q_{\mathcal{A}})$$

$${\eta }_{\mathrm{mat}}(q_{\mathcal{A}})=$$

$$\alpha =$$

$$+2qn+2q,$$

$$\beta =$$

$$\frac{{\eta }_{\mathrm{inv}}(q_{\mathcal{A}})}{{\eta }_{\mathrm{mat}}(q_{\mathcal{A}})}\le \frac{18}{79}$$

$$(2q_{\mathcal{A}}^2+12q_{\mathcal{A}}+\frac{10}{3})(\alpha q_{\mathcal{A}}+\beta )-(\frac{2}{3}{q}_{\mathcal{A}}^2+6q_{\mathcal{A}}+\frac{10}{3})\alpha q_{\mathcal{A}}$$

$$(\frac{4}{3}{q}_{\mathcal{A}}^2+6q_{\mathcal{A}})\alpha q_{\mathcal{A}}+(2q_{\mathcal{A}}^2+12q_{\mathcal{A}}+\frac{10}{3})\beta >0.$$

$$N=2,n=1,m\ge 1,$$

$$N=2,n\ge 2,m\ge 1,$$

$$N\ge 3,n\ge 1,m\ge 1,$$

$$1328m^2+3888mn+1296n^2-5256m+486n-1742\ge 0.$$

Full text

Introduction

Model predictive control (MPC) is an established method for linear, multivariate systems with constraints. MPC is computationally expensive, however. MPC controllers are often too complex for an implementation on standard industrial hardware such as programmable logic controllers (PLC).

Various approaches to reducing the computational effort of MPC have been investigated. In [2] and [24], parametric programming was used to solve the optimization problem offline. Code generation tools ([28, 12, 19]) produce a problem-tailored implementation of the MPC problem. In [21] and [15], methods for an online simplification of the optimal control problem have been presented. Event-based approaches, other than the one discussed here, have been proposed in [7, 27, 16, 14, 3] and the references therein (see also the introduction in [11]).

Tailored MPC implementations for hardware with limited resources such as PLCs have also been studied. A suboptimal parametric programming approach has been used in [26]. The degree of suboptimality was chosen so that the resulting controller satisfied resource restrictions of the PLC. Parametric programming was also used in [22], where a binary search tree is instrumental to an implementation with limited memory. In [13], the use of various quadratic program (QP) solvers for online MPC on PLC was investigated. The results were illustrated with a temperature control problem. In [17] primal-dual first-order methods were used to speed up solving the optimization problem and the proposed method was tested with a hardware-in-the-loop (HIL) simulation.

We use an event-triggered MPC scheme [5] that we claim to be particularly suitable for networked MPC. The central idea is as follows: Solving a linear-quadratic optimal control problem (OCP) at the current state is usually assumed to provide the optimal control action for this particular state. It is well-known, however, that this solution additionally provides an affine state feedback law that is optimal on a state space polytope [2, 24]. As long as the system state remains in this polytope, the feedback law can be evaluated at very low computational effort. Consequently, lean local nodes such as PLCs, directly attached to controlled systems, suffice to evaluate affine feedback laws as long as the system stays in the polytope of validity. Upon leaving the polytope, a powerful central node is called to solve an OCP and to provide a new feedback law and polytope. This setup obviously requires affine feedback laws (and linear inequalities that define polytopes) to be sent across the network. Since bandwidth and computational resources are limited, it is crucial to investigate alternatives to sending the feedback law and polytope (i.e., real matrices and vectors) across the network. The present paper extends earlier theoretical results [5] by new results on the computational and bandwidth requirements of the alternatives (Propositions 4-6), by details on the implementation (Section 4), and by a validation with an actual implementation on industrial hardware (Section 5).

We introduce the event-triggered MPC scheme in Section 2. Section 3 discusses the computational effort and bandwidth requirements of the event-triggered controller. Section 4 addresses the implementation on industrial hardware. Some results obtained with this implementation, which corroborate our theoretical complexity analysis, are discussed in Section 5. A conclusion is drawn in Section 6.

Notation

Let $M_{\mathcal{A}}$ refer to the submatrix of the matrix $M$ with the row indices listed in the index set $\mathcal{A}$, where all index sets are assumed to be ordered. A polytope $\mathcal{P}$ is defined to be the intersection of a finite number of halfspaces, i.e., for any polytope $\mathcal{P}$ there exist $T\in\mathbb{R}^{q\times n},\ d\in\mathbb{R}^{q}$ such that

[TABLE]

Event-triggered networked MPC

Problem class

We consider linear discrete-time MPC problems, where the optimal control problem

[TABLE]

with state variables $x\in\mathbb{R}^{n}$ and input variables $u\in\mathbb{R}^{m}$, horizon $N\in\mathbb{N}>1$, weighting matrices $P\in\mathbb{R}^{n\times n}$, $P\succ 0$, $Q\in\mathbb{R}^{n\times n}$, $Q\succeq 0$ and $R\in\mathbb{R}^{m\times m}$, $R\succ 0$ is solved. The system $(A,B)$ with matrices $A\in\mathbb{R}^{n\times n}$ and $B\in\mathbb{R}^{n\times m}$ is assumed to be stabilizable. States and inputs are subject to constraints (2c), (2d), (2e) with compact polytopes $\mathcal{U}\subset\mathbb{R}^{m}$, $\mathcal{X}\subset\mathbb{R}^{n}$ and $\mathcal{T}\subseteq\mathcal{X}$ that contain the origin in their interiors.

The MPC problem (2) can be stated as a quadratic program (QP) of the form [18, Chap. 3]:

[TABLE]

with $U=(u(0)^{\prime},\dots,u(N-1)^{\prime})^{\prime}\in\mathbb{R}^{mN}$, $F\in\mathbb{R}^{n\times mN}$, $H\in\mathbb{R}^{mN\times mN}$, $G\in\mathbb{R}^{q\times mN}$, $w\in\mathbb{R}^{q}$, $E\in\mathbb{R}^{q\times n}$, where $\mathcal{Q}=\left\{1,\ldots,q\right\}$ and $q$ is the number of constraints in (2) and (3). $H$ is positive definite since $Q\succeq 0$ and $R\succ 0$ [18, p. 76]. The set of all feasible states is denoted by $\mathcal{X}_{f}$ and assumed to be nonempty. Since $H\succ 0$, there exists a unique optimal solution for every $x\in\mathcal{X}_{f}$ that we denote $U^{\star}(x)$. There exist polytopes $\mathcal{P}_{1},...,\mathcal{P}_{p}$ and affine control laws $K_{i}x+b_{i}$, $K_{i}\in\mathbb{R}^{mN\times n}$, $b_{i}\in\mathbb{R}^{mN}$, $i=1,...,p$, such that $\cup_{i}^{p}\mathcal{P}_{i}=\mathcal{X}_{f}$ and

[TABLE]

and $U^{\star}(x)$ is a continuous function on $\mathcal{X}_{f}$ [2, 24]. We stress that we exploit the structure of the explicit solution (4) without ever calculating it entirely.

Event-triggered controller

Equation (4) suggests the simple networked MPC setup that was already briefly described in the introduction: Solving (2) or (3) for the current state $x(k)\in\mathbb{R}^{n}$ is usually understood to result in the optimal input signal for this particular point $x(k)$ in the state space. It follows from (4), however, that the solution to (3) for a point $x\in\mathcal{X}_{f}$ not only provides the optimal input signal at this point, but an affine law $K_{i}x+b_{i}$ and its polytope of validity $\mathcal{P}_{i}$ (see Lemma 1 below). For all $x\in\mathcal{P}_{i}$, this affine law yields the same optimal input as solving the optimal control problem (2) or (3). Since the computational effort for evaluating $K_{i}x+b_{i}$ and checking whether $x\in\mathcal{P}_{i}$ is much smaller than that for solving (2) or (3), the optimal solution can be determined on a lean local node such as a PLC with $K_{i}x+b_{i}$ and $\mathcal{P}_{i}$. A computationally powerful central node (the industrial PC in the particular setup investigated here) is only required to solve the OCP whenever the current polytope is left.

The calculations necessary to obtain the affine law $K_{i}x+b_{i}$ and its polytope of validity $\mathcal{P}_{i}$ from $U^{\star}(x)$ are summarized in the following Lemma for ease of reference.

Lemma 1

Let $U^{\star}(x)$ be the solution to (3) for an arbitrary $x\in\mathcal{X}_{f}$, and let $\mathcal{A}$ and $\mathcal{I}$ refer to the sets of active and inactive constraints

[TABLE]

Let $\mathcal{M}=\{1,\dots,m\}$ and assume the matrix $G_{\mathcal{A}}$ has full row rank. Then

[TABLE]

with $S=E+GH^{-1}F^{\prime}$, $S\in\mathbb{R}^{q\times n}$ define the optimal control law and the polytope

[TABLE]

A concise proof of Lemma 1, which is an immediate corollary to the results in [2], is given in [16].

For later use we introduce $q_{\mathcal{A}}=|\mathcal{A}|$ and

[TABLE]

The function $e(x)$ is used to check whether the control law (7a) can be reused to compute the control input for the current state $x(k)$ ($e(x(k))=0$), or a new optimal control law has to be calculated ($e(x(k))=1$). We introduce the term trigger event, or in short event to describe any instant for which $e(x(k))=1$. The inequalities in (7b) are validated row by row in the implementation. If all inequalities hold, no event occurred (i.e., $e(k)=0$).

Event-triggered networked MPC

Consider the simple star-shaped network architecture depicted in Fig. 1, where multiple local control nodes are connected to a single central node via network connections. Throughout this paper we call this structure an event-triggered networked model predictive controller in contrast to classical MPC, where the OCP is solved and the control input is applied to the system in every time step.

In case of an event, a control law and polytope must be sent across the network. The simplest way of doing so is to send matrices $K_{\mathcal{M}}^{\star},b_{\mathcal{M}}^{\star},T^{\star}$ and $d^{\star}$, that represent the control law (7a) and polytope (7b). However, equations (5-7b) suggest to distribute computations between local and central node such that the amount of transmitted data may be reduced. It is our aim to identify the approach that results in the best trade-off between transmitted data and computations on the local node, and to propose an implementation scheme for this particular variant.

Computational- and network load in event-triggered networked MPC

Measuring transmitted data and computational cost

The remainder of the paper is based on the following assumptions and hypotheses H1-H3: (H1) Bandwidth requirements can be determined by counting the number of bits transferred per event, where real numbers are transmitted using IEEE 754 half precision [1] ($\lambda=16\frac{\text{bit}}{\text{real}}$ below). (H2) Computational effort can be measured by counting floating point operations, where multiplications and summations cost one floating point operation (flop), and divisions cost ten flops [25]; the cost of some recurring operations is stated in Table 1 [10, chap. 1.1.15] [8, p.12]. (H3) The central node is assumed to carry out all computations instantaneously.

Let the term computational load on the local node refer to the number of flops necessary to determine the control law and polytope (7). Since all variants discussed below have the evaluation of the control law (7a) and condition (8) in common, we neglect the flops required for (7a) and (8).

Trade-off between transmitted data and local calculations

A control law and polytope must be sent across the network whenever an event (8) occurs. It is evident from (7) that the matrices $K_{\mathcal{M}}^{\star},b_{\mathcal{M}}^{\star},T^{\star}$ and $d^{\star}$ define the control law and polytope. It suffices to send much less data, however. This is stated more precisely in Lemma 2, which was proved in [5].

Lemma 2

Let $x\in\mathcal{X}_{f}$ be arbitrary and assume the optimal control problem (3) has been solved for this initial condition. Then the feedback law (7a) and its polytope of validity (7b) are determined by any of the following:

[TABLE]

Equation (9b) may tempt the reader to think we use the input signals $u^{\star}(0),\dots,u^{\star}(N-1)$ of $U^{\star}$, which is open-loop but in general not closed-loop optimal. We stress again this is not the case, but we use the feedback law (7a) and polytope (7b), because this feedback law is closed-loop optimal on the polytope. Because we determine a new law and polytope whenever the current polytope is left, the same closed-loop optimal sequence as solving (2) in every time step results. Therefore, the same closed-loop optimal system behavior is obtained as solving (2) or (3) in every time step [16]. The claim about (9b) in Lemma 2 holds, because $U^{\star}$ defines $\mathcal{A}(x)$ according to (5) and $\mathcal{A}$ defines $K^{\star}_{\mathcal{M}}$, $b^{\star}_{\mathcal{M}}$, $T\star$, $d^{\star}$ and thus the law and polytope according to (6).

Lemma 2 implies the following alternative approaches A1-A4 for transmitting the optimal control law (7a) and its polytope (7b).

• A1

Send $q$ bits

[TABLE]

$i=1,\dots,q$, from the central to the local node, which is equivalent to sending $\mathcal{A}$ and $\mathcal{I}$ according to (5). Evaluate (6) on the local node to determine $K_{\mathcal{M}}^{\star}$, $b_{\mathcal{M}}^{\star}$, $T^{\star}$ and $d^{\star}$.

• A3

Determine

[TABLE]

on the central node and send it in addition to and redundantly to (10). Evaluate  (6) to determine $K_{\mathcal{M}}^{\star}$, $b_{\mathcal{M}}^{\star}$, $T^{\star}$ and $d^{\star}$ as in A1 on the local node. Use the inverse $\Phi$ instead of determining it on the local node.

• A5

Send $U^{\star}$ from the central to the local node. Evaluate (5) to determine $\mathcal{A}$ and $\mathcal{I}$, then evaluate (6) to determine $K_{\mathcal{M}}^{\star}$, $b_{\mathcal{M}}^{\star}$, $T^{\star}$ and $d^{\star}$ on the local node.

• A7

Determine $K_{\mathcal{M}}^{\star}$, $b_{\mathcal{M}}^{\star}$, $T^{\star}$, $d^{\star}$ on the central node and send them to the local node.

These four approaches obviously differ with respect to the amount of data sent across the network and the computational effort on the local node. More precisely, we have the following results, which are taken from [5].

Lemma 3

Assume H1-H4 hold. Furthermore, assume (2) or equivalently (3) has been solved for an $x\in\mathcal{X}_{f}$ on the central node. Then the number of bits that need to be transmitted in approaches A1-A4 and the number of flops required on the local node for approaches A1-A4 are as follows:

[TABLE]

Based on the technical results stated in Lemma 3, it is now easy to compare approaches A1-A4. We compare the amount of transmitted data first.

Proposition 4

Let all assumptions be as in Lemma 2. The following statements hold for the number of transmitted bits in approaches A1-A7:

• (i)

Approach A1 requires as much as, or less data than approach A3 to be sent across the network.

If we additionally assume all constraints (2c, 2d, 2e) to be box constraints (i.e., if $q=2mN+2nN+2n$), the following statements hold:

• (ii)

Approaches A1-A5 never require more data to be transmitted than approach A7.

• (iii)

Approach A5 requires more data to be sent across the network than approach A1, if

[TABLE]

where $\lambda$ is the number of bits per real number.

• (iv)

Approach A5 requires more data to be sent across the network than approach A3, if

[TABLE]

Proof: Statement (i) is obvious, since A1 requires to transmit $\gamma$ and approach A3 requires to transmit $\gamma$ and $\Phi$. Statement (ii) can be shown by substituting $q=2mN+2nN+2n$ (box constraints) into the number of transmitted bits stated in Lemma 3 and subtracting the number of bits for A1-A5 from those for A7. Since the resulting expressions are always positive, statement (ii) holds. In order to prove (iii) we need to show that (12) implies

[TABLE]

For $q=2mN+2nN+2n$ (box constraints), this is equivalent to

[TABLE]

It is easy to show $N>1$ implies $2\frac{n}{m}\left(1+\frac{1}{N}\right)\leq 3\frac{n}{m}$. Therefore,

[TABLE]

implies (15), i.e., (14) for box constraints. Rearranging (16) yields the desired condition (12). Part (iv) can be shown in the same way as (iii). $\hfill\blacksquare$

Since states and inputs are constrained by an upper and lower bound in many technical systems, the introduction of box constraints in part (ii) - (iv) of Proposition 4 is not a severe restriction.

Proposition 4 reveals the following partial order of the approaches with respect to the number of transmitted bits: Approach A1 requires the smallest number of transmitted bits, followed by approaches A3 and A7. Approach A5 results in less transmitted data than A1 if inequality (12) is not fulfilled and less transmitted data than approach A3 if inequality (13) is not fulfilled. It remains to compare the computational effort on the local node based on Lemma 2.

Proposition 5

Let all assumptions be as in Lemma 2. The following statements hold for the computational effort on the local node.

• (i)

The number of flops in approach A3 never exceeds those needed in approach A1.

• (ii)

Approach A5 always requires more computations than approach A1.

• (iii)

Approach A7 requires no operations on the local node to determine the matrices for the control law and polytope.

Proof: Statement (iii) is obvious from the zero in the last row and column of the table in Lemma 3. Statements (i) and (ii) can be shown by subtracting the number of flops given in the third column of the table in Lemma 3 for the respective pair of approaches. Specifically, this difference reads $\frac{2}{3}q_{\mathcal{A}}^{3}+q_{\mathcal{A}}^{2}(2mN+5)+\frac{10}{3}q_{\mathcal{A}}$ for approach A1. Since this expression is positive, approach A1 never requires less operations than approach A3. Statement (ii) can be shown in the same way. $\hfill\blacksquare$

None of the approaches A1-A7 results in both, the smallest number of transmitted bits and the smallest computational effort. Since the number of transmitted bits in approach A3 depends on $q_{\mathcal{A}}$ (see (13)), and since approach A5 depends on (12), the best trade-off results for transmitting the active set of constraints in approach A1. We examine A1 in more detail in the remainder of the paper.

Computational effort of approach A1

While only $q_{\mathcal{A}}$ bits need to be transmitted in approach A1, additional computations must be carried out on the local node. Most importantly, the inverse of $G_{\mathcal{A}}H^{-1}(G_{\mathcal{A}})^{\prime}\in\mathbb{R}^{q_{\mathcal{A}}\times q_{\mathcal{A}}}$ has to be determined on the local node, which requires on the order of $q_{\mathcal{A}}^{3}$ operations (see, e.g., [8], p.12). However, the absolute number of operations for this inversion always remains small compared to the number of operations necessary for evaluating (6), even though the latter scales linearly only in $q_{\mathcal{A}}$. This statement is made more precise in the following proposition for the special case of box constraints. Note that the proposition is conservative in the sense that the computational effort for the inversion is even smaller in practice (see Section 4).

Proposition 6

Let $\eta_{\rm inv}(q_{\mathcal{A}})$ and $\eta_{\rm mat}(q_{\mathcal{A}})$ refer to the number of flops necessary to carry out the matrix inversion (11) and and all other matrix operations in (6), respectively. Assume the constraints (2c)-(2e) to be box constraints, which implies $q=2mN+2nN+2n$, and assume $N>1$. Then

[TABLE]

Proof: The inversion in (6) requires

[TABLE]

flops according to [8, p.12]. Subtracting this figure from the overall number of operations for A1 as stated in Lemma 3 yields

[TABLE]

which equals the number of flops for all matrix operations but the inversion in A1. It is easy to show that $m>0$, $N>1$, $q>0$ and $n>0$ imply $\alpha>0$ and $\beta>0$. We now show that if

[TABLE]

holds for a $\bar{q}_{\mathcal{A}}$, it holds for all $q_{\mathcal{A}}<\bar{q}_{\mathcal{A}}$. This can be seen by extending $q_{\mathcal{A}}\in\mathbb{N}$ to a real variable $q_{\mathcal{A}}\in\mathbb{R}$ and showing $\eta_{\rm inv}/\eta_{\rm mat}$ is strictly increasing by proving $\left.\frac{d}{dq_{\mathcal{A}}}\!\left(\frac{\eta_{\rm inv}(q_{\mathcal{A}})}{\eta_{\rm mat}(q_{\mathcal{A}})}\right)\!>0\right.$. Applying the quotient rule yields a fraction with numerator

[TABLE]

and denominator $\eta_{\rm mat}^{2}(q_{\mathcal{A}})$. Since the denominator is positive, it suffices to show (19) is positive, which can equivalently be stated as

[TABLE]

This condition obviously holds for $\alpha>0$, $\beta>0$ and $q_{\mathcal{A}}>0$. It remains to be shown that (18) is positive for the largest possible number of concurrently active constraints $\bar{q}_{\mathcal{A}}$. This number reads $q_{\mathcal{A}}=\bar{q}_{\mathcal{A}}=mN$ for box constraints with $q=2mN+2nN+2n$. All combinations of $N>1$, $n>0$ and $m>0$ are covered by the three cases

[TABLE]

where (i) and (ii) together cover $N=2$, $n>0$, $m>0$. Case (i), i.e., substituting $N=2$ and $n=1$ into (18) and rearranging yields the condition $1328m^{2}-1368m+40\geq 0$, which obviously holds for all $m\geq 1$. Case (ii), i.e., substituting $N=2$ into (18) and rearranging yields the condition

[TABLE]

Now $n\geq 2$ implies $3888mn-5256m\geq 0$ and $1296n^{2}+486n-1742\geq 0$ for all $m$. Adding these two inequalities to $1328m^{2}\geq 0$, which also holds for all $m$, shows that condition (20) holds. Case (iii) can be shown by noting that (18) can equivalently be stated as

[TABLE]

for arbitrary $N$, $m$ and $n$. It is easy to see that $N\geq 3$ and $m\geq 1$ imply the first three terms on the left hand side of (21) are larger than the first term on the right hand side for all $\left.n\geq 1\right.$. Similarly, the 4th and 5th term on the left hand side of (21) together dominate the second term on the right hand side, and the last term on the left hand side dominates the remaining two terms on the right hand side. Since the remaining terms on the left hand side are positive, condition (21) holds. $\hfill\blacksquare$

Implementation on industrial hardware

We use a SIMATIC IPC847C industrial computer (IPC for short) with an Intel Core i7-610E CPU and 8GB memory111See http://www.siemens.com for more information about the IPC, PLC and I/O modules (checked on October 7, 2016). as the central node. The IPC runs a QNX Neutrino 6.5.0 real-time operating system and is programmed in ISO C++. The local node consists of a Siemens S7-400 PLC equipped with 30MB memory, an SM 331 module, and an SM 332 analog I/O module22footnotemark: 2. We use statement lists (STL) to program the PLC. The system that is controlled by the PLC is simulated on a dSpace222 See https://www.dspace.com for more information (checked on October 7, 2016). HIL system. System states and inputs are exchanged via analog I/Os. We refer to this equipment as HIL setup for brevity.

Computation times are measured on the PLC by setting a digital output to high every time the computational task starts, and to low when it finishes. Closer inspection of these measurements reveals that the digital I/O modules of the PLC have a small random delay. We repeat measurements ten times to reduce the influence of this delay. This results in standard deviations of less than 3% for all reported data.

The communication between the nodes is performed by a PROFIBUS network, but we stress that any other bus system can be used instead. The network is only used when the local node triggers the central node to calculate a new affine control law and polytope.

Implementation details

We use a simple code generator to generate the required STL code from the problem (2) in Matlab. The code generator precomputes matrices $H^{-1},G,S,F^{\prime},w$ offline. The evaluation of (6) exploits common subexpressions. The inversion of $G_{\mathcal{A}}H^{-1}(G_{\mathcal{A}})^{\prime}$ is replaced by a LU decomposition with pivoting and forward-backward substitutions. The code generator also generates a C++ header file for the IPC which provides matrices of the QP (3). We use OOQP (see [9]) to solve QPs on the IPC.

Since data transmission and calculations (3), (6) may be necessary to obtain the control input $u$, a delay between measuring the state and applying the control input may occur. However, the delay only occurs if $e(k)=1$ in (8). Consequently, depending on whether $e(k)=0$ or $e(k)=1$ in (8), the computation of the control input differs. We explain these two cases in more detail. Let $x(k)$ and $x(k+1)$ refer to the current and subsequent state of the closed-loop system, and let $\mathcal{P}^{\star}_{k}$ be a polytope with $x(k)\in\mathcal{P}^{\star}_{k}$. While $x(k)$ is known from a measurement, $x(k+1)$ is predicted at sampling interval $k$ with the system model (2b). Either one of the following two cases applies in every time step:

(i) If $x(k+1)\in\mathcal{P}^{\star}_{k}$ , the control law and polytope in the subsequent time step $k+1$ remain the same as in time step $k$. The control input $u(k+1)$ is calculated at the beginning of sampling interval $k+1$ on the local node and is applied to the system. In this case no communication with the central node is necessary.

(ii) If $x(k+1)\not\in\mathcal{P}^{\star}_{k}$, the QP (3) must be solved on the central node and a new control law and polytope must be calculated locally. The time between two consecutive sampling intervals $k$ and $k+1$ can be used to compute the control law and polytope (7) that are required in time step $k+1$: We transfer the predicted state $x(k+1)$ to the central node, solve the quadratic program (3) and transfer data according to approaches A1-A7 back to the local node. Once the information is available on the local node, matrices of the control law and polytope are calculated. By solving (3) for the *predicted state * $x(k+1)$, the control law and polytope are obtained for time step $k+1$. The control input $u(k+1)$ is calculated at the beginning of the next sampling interval by evaluating the predicted control law.

The flow chart in Fig. 2 visualizes computations and communication in the HIL setup. In case (i) all calculations are done on the local node (yellow box). In case (ii), the central node (blue box) solves the QP (3).

Validation of the implementation

We validated the implementation on the PLC and the IPC by regulating the simple double integrator from [24] to the origin for 100 random initial states and comparing the state space trajectories that resulted from the HIL setup to reference trajectories simulated in Matlab. A mean squared error $\frac{1}{M}\sum\nolimits_{i=1}^{M}\left\lVert x(i)_{\rm Matlab}-x(i)_{\rm HIL}\right\rVert=0.0071$ resulted, where a total of $M=4051$ time steps arose for the 100 random initial conditions. The error results from analog to digital and digital to analog conversion errors, and random noise on the analog signal.

Results

We measure computation times of approach A1 running on the HIL setup described in Section 4. The results corroborate Proposition 6.

Example system

Consider the mechanical system with four oscillating masses connected by springs sketched in Fig. 3. Three inputs ($u_{l},l=1,2,3$) are used to control tension between the masses. Masses and spring constants have the value one and there is no damping. The resulting continuous time system

[TABLE]

is discretized with zero order hold and a sampling time of $T_{s}=0.5$ seconds. State variables $x\in\mathbb{R}^{8}$ and input variables $u\in\mathbb{R}^{3}$ are subject to the constraints $x_{i}\in\left[-4,4\right],i=1,...,8$ and $u_{l}\in\left[-0.5,0.5\right],l=1,2,3$. Cost function matrices are chosen to be $Q=I^{n\times n}$ and $R=I^{m\times m}$, and $P$ is set to the solution of the discrete-time algebraic Riccati equation. The horizon of the MPC problem is chosen to be $N=10$. The number of simultaneously active constraints is bounded above according to

[TABLE]

Hardware in the loop test

Approach A1 provides the best trade-off of data transmission to computational cost on the local node. We conduct measurements of the computation times on the local node with the HIL setup described in Section 4 to analyze approach A1 in more detail. Specifically, we report actual computational times of the HIL system for the following tasks to corroborate Proposition 6:

• a)

transferring the current system state to the central node, solving the QP on the central node, and sending $\gamma$ as defined in (10) to the local node,

• b)

evaluation of (6) on the local node,

• c)

computing the control input and evaluating condition (7) on the local node.

We assume that flops and computation times are related by a constant factor, i.e.

[TABLE]

for some $\alpha\in\mathbb{R}^{+}$.

The computational times that result for 100 random feasible initial conditions are summarized in Fig. 5. The times required for tasks a, b and c can be approximated by a constant, a cubic function, and a constant, respectively. The time required for task c is obviously small compared to those of tasks a and b, which demonstrates that the control input can be computed with very low computational effort if $e(x)=0$. The time required for task b dominates the computational time required on the local node for all but very small $q_{\mathcal{A}}$.

As suggested by Proposition 6, we further analyze the computation time $\Delta t_{b}$ for task b by splitting it into the time $\Delta t_{{\rm inv}}$ necessary to calculate the inverse of $G_{\mathcal{A}}H^{-1}(G_{\mathcal{A}})^{\prime}$ and the time $\Delta t_{{\rm mat}}$ required for all remaining matrix operations in (6). Figure 5 shows the two contributions

[TABLE]

Substituting

[TABLE]

which holds according to Proposition 6 and with (24), results in

[TABLE]

This can be understood as an upper bound on $\Delta t_{b}$, which is also shown in Fig. 5. It is evident from Fig. 5 that the absolute contribution of $\Delta t_{{\rm inv}}$ remains small despite the cubic dependence of $\Delta t_{{\rm inv}}$ on $q_{\mathcal{A}}^{3}$. We conclude approach A1 is the best one of the four variants, because the inversion of $G_{\mathcal{A}}H^{-1}(G_{\mathcal{A}})^{\prime}$ required in A1 but not in A3 results in only a small additional computational effort on the local node. Note that this result does not only apply to the example, but the bound (17) established in Proposition 6 applies in general.

Conclusion and outlook

We analyzed the practical usefulness of a networked event-triggered MPC approach by implementing it on industrial hardware. The proposed approach combines a powerful compute node, which solves optimal control problems on demand, with lean local nodes that control local systems with simple affine, regional optimal feedback laws. The local nodes use their feedback laws as long as they are optimal and otherwise request a new law from the central node. Several options for an implementation exist that differ with respect to their balance of local computational effort and network bandwidth requirements. We extended earlier theoretical results to determine the best variant among these choices and corroborated the theoretically expected features of the selected variant with an implementation on industrial control hardware.

An extension to tube-based robust MPC is possible (see [23] for a simulation-based analysis). We believe an extension in this direction is more of theoretical interest, because the case treated in the present paper inherits the robustness of the regional approach [5]. This aspect is currently under investigation.

The proposed approach is based on the idea that once the set of active constraints is known for the current system state, an optimal feedback law and the polytope on which it is optimal are known (see Lemma 2). This also holds in nonlinear MPC with discrete-time system models [20, 6]. It therefore suffices to send the active sets across the network in nonlinear MPC. However, the local nodes must solve nonlinear systems of equations in this case, while affine control laws only need to be evaluated in the linear MPC case. It remains to investigate how the computational effort on the local nodes can be reduced, for example by sending redundant information as in approach A3.

Finally, future work must combine the proposed approach with methods for predicting the sequence of affine laws along the closed-loop trajectory [4]. This will help reducing the number of optimal control problems further in particular far from the origin, where the affine laws are not reused frequently.

Acknowledgment

Support by the Deutsche Forschungsgemeinschaft (DFG) under the grant MO 1086/15-1 is gratefully acknowledged.