Control Synthesis for Permutation-Symmetric High-Dimensional Systems With Counting Constraints

TL;DR

This paper introduces a method for synthesizing correct controllers for high-dimensional, symmetric systems with counting constraints by leveraging symmetry and aggregate abstractions, enabling scalability to tens of thousands of states.

Contribution

It presents a novel approach combining aggregate abstraction and linear inequality formulation to handle large, symmetric systems with complex counting constraints.

Findings

Successfully synthesizes controllers for systems with tens of thousands of states.

Demonstrates scalability and correctness of the approach.

Effectively manages symmetry to reduce computational complexity.

Abstract

General purpose correct-by-construction synthesis methods are limited to systems with low dimensionality or simple specifications. In this work we consider highly symmetrical counting problems and exploit the symmetry to synthesize provably correct controllers for systems with tens of thousands of states. The key ingredients of the solution are an aggregate abstraction procedure for mildly heterogeneous systems and a formulation of counting constraints as linear inequalities.

Figures

Tables

Table 1. TABLE I: Given feasibility or infeasiblity of ( 24 ) in various configurations, this table lists the inferences that can be made according to results in Section VI .

LP/ILP	Feasible?	Cycle set	$T$	Result	Why?
ILP	Yes	Any	Any	Solution	Theorem 3
LP	Yes	Any	Any	Approximate solution	Theorem 7,8
ILP	No	$\{C:\left|C\right|\le \left|𝒬\right|\left(\frac{\left|𝒬\right|+N-1}{N}\right)\}$	$\left(\frac{\left|𝒬\right|+N-1}{N}\right)$	No solution exists	Theorem 5
LP	No with $ϵ$-relaxation	$\{C:C\text{ simple}\}$	$\frac{(\mathrm{diam}{\left(G\right)}^2+1)N}{ϵ}$	No solution exists	Theorem 6

Table 2. TABLE II: Average solution times over 10 randomized trials

N	${10}^2$	${10}^3$	${10}^4$	${10}^5$	${10}^6$	${10}^7$	${10}^8$	${10}^9$
Time (s)	7.8	11.1	12.0	12.8	11.5	13.2	13.0	11.0

Table 3. TABLE III: Parameter values for the two classes of subsystems. Note that the time discretizations need to be identical for concurrent execution.

Parameter	Class 1	Class 2
Nominal $[a,b,{\theta }_a,P_m]$	$[2,2,32,5.6]$	$[2.2,2.2,32,5.9$]
Space discretization $\eta$	0.002	0.0015
Time discretization $\tau$	0.05	0.05
Error bound $\overline{\delta }$	0.025	0.025

Equations

$$\frac{\mathrm{d}}{\mathrm{d}t}{\mathbf{x}}_n(t)=f_{{\sigma }_n(t)}\left({\mathbf{x}}_n(t),{\mathbf{d}}_n(t)\right),{\sigma }_n(t)\in \{\text{on},\text{off}\},$$

$$\text{xn(t)∈[a,a] for all n∈[N] and all t∈R+},$$

$$\begin{array}{rl}& \text{For all t∈R+, at least K and at most K}\\ & \text{ of the TCLs are in mode on.}\end{array}$$

$$\sum _{n\in [N]}\text{mathds}{1}_X\left(q_n,{\mu }_n\right)\le R.$$

$$\sum _{n\in [N]}\text{mathds}{1}_{\mathbb{R}\times \{\text{on}\}}\left({\mathbf{x}}_n(t),{\sigma }_n(t)\right)\le \overline{K},$$

$$\sum _{n\in [N]}\text{mathds}{1}_{\mathbb{R}\times \{\text{off}\}}\left({\mathbf{x}}_n(t),{\sigma }_n(t)\right)\le N-\underline{K}.$$

$$\sum _{n\in [N]}\text{mathds}{1}_{[\underline{a},\overline{a}{]}^C\times \{\text{on},\text{off}\}}\left({\mathbf{x}}_n(t),{\sigma }_n(t)\right)\le 0.$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{\mathbf{x}}_i(t)=f_{{\sigma }_i(t)}({\mathbf{x}}_i(t),{\mathbf{d}}_i(t)),{\sigma }_i:\mathbb{R}\to \mathcal{U},$$

$$\Vert {\varphi }_{\mu }(t,x,0)-{\varphi }_{\mu }(t,y,0){\Vert }_{\infty }\le {\beta }_{\mu }\left(\Vert x-y{\Vert }_{\infty },t\right).$$

$$\Vert f_{\mu }(x,d)-f_{\mu }(x,0){\Vert }_{\infty }\le {\bar{\delta }}_{\mu },\forall d\in \mathcal{D}.$$

$$\sum _{n\in [N]}\text{mathds}{1}_{X_l}\left({\mathbf{x}}_n(s),{\sigma }_n(s)\right)\le R_l,\forall s\in \mathbb{N},\forall l\in [L].$$

$${\Sigma }_{\tau }=(\mathcal{X},\mathcal{U},\underset{\tau }{⟶},{\text{Id}}_{{\mathbb{R}}^d}),$$

$$(N,{\Sigma }_{\tau },\{{\mathbf{x}}_n(0){\}}_{n\in [N]},\{X_l,R_l{\}}_{l\in [L]})$$

$${\Sigma }_{\tau ,\eta }=\left({\kappa }_{\eta }(\mathcal{X}),\mathcal{U},\underset{\tau ,\eta }{⟶},{\text{Id}}_{{\mathbb{R}}^d}\right),$$

$${\beta }_{\mu }(ϵ,\tau )+\frac{{\overline{\delta }}_{\mu }}{K_{\mu }}\left(e^{K_{\mu }\tau }-1\right)+\frac{\eta }{2}\le ϵ,$$

$${\mathcal{G}}^{+ϵ}(X^{\mathcal{X}}\times X^{\mathcal{U}})={\kappa }_{\eta }\left(X^{\mathcal{X}}\oplus {\mathcal{B}}_{\infty }(0,ϵ)\right)\times X^{\mathcal{U}},$$

$${\mathcal{G}}^{-ϵ}(X^{\mathcal{X}}\times X^{\mathcal{U}})={\kappa }_{\eta }\left(X^{\mathcal{X}}\ominus {\mathcal{B}}_{\infty }(0,ϵ)\right)\times X^{\mathcal{U}}.$$

$$\left(N,{\Sigma }_{\tau ,\eta },\{{\kappa }_{\eta }({\mathbf{x}}_n(0)){\}}_{n\in [N]},\{{\mathcal{G}}^{+ϵ}(X_l),R_l{\}}_{l\in [L]}\right)$$

$$(N,{\Sigma }_{\tau },\{{\mathbf{x}}_n(0){\}}_{n\in [N]},\{X_l,R_l{\}}_{l\in [L]}).$$

$$\left(N,{\Sigma }_{\tau ,\eta },\{{\kappa }_{\eta }({\mathbf{x}}_n(0)){\}}_{n\in [N]},{\left\{{\mathcal{G}}^{-(ϵ+\frac{\eta }{2})}(X_l),R_l\right\}}_{l\in [L]}\right)$$

$$\left(N,{\Sigma }_{\tau },\{{\mathbf{x}}_n(0){\}}_{n\in [N]},\{X_l,R_l{\}}_{l\in [L]}\right).$$

$$\sum _{n\in [N]}\text{mathds}{1}_{[\underline{a}+ϵ,\overline{a}-ϵ{]}^C\times \{\text{on},\text{off}\}}\left({\xi }_n(s),{\sigma }_n(s)\right)\le 0,$$

$$w_q(s+1)=\sum _{\mu \in \mathcal{U}}\sum _{q^{\prime }\in {\mathcal{N}}_q^{\mu }}{r}_{q^{\prime }}^{\mu }(s),q\in \mathcal{Q},$$

$$r_q^{\mu }(s)\ge 0,\sum _{\mu \in \mathcal{U}}{r}_q^{\mu }(s)=w_q(s),$$

$$\Gamma :\mathbf{w}(s+1)=B\mathbf{r}(s),$$

$$\mathcal{W}$$

$$\mathcal{R}(\mathbf{w})$$

$$C=(q_0,{\mu }_0,q_1)(q_1,{\mu }_1,q_2)\dots (q_{I-1},{\mu }_{I-1},q_0),$$

$${\alpha }^{\circlearrowleft s}(i)=\alpha \left((i-s)\mathrm{mod}|C|\right).$$

$${⟨C,{\alpha }^{\circlearrowleft s}⟩}^X=\sum _{i\in [|C|]}\{\begin{array}{ll}{\alpha }^{\circlearrowleft s}(i)& \text{if}(q_i,{\mu }_i)\in X,\\ 0& \text{otherwise}.\end{array}$$

$${\mathrm{maxcnt}}^X\left(C,\alpha \right)=\underset{s\in \mathbb{N}}{\max }{⟨C,{\alpha }^{\circlearrowleft s}⟩}^X.$$

$$\Vert B_C^X\alpha {\Vert }_{\infty }={\Vert \left[\begin{array}{ccccc}0& 1& 1& 1& 0\\ 1& 1& 1& 0& 0\\ 1& 1& 0& 0& 1\\ 1& 0& 0& 1& 1\\ 0& 0& 1& 1& 1\end{array}\right]\left[\begin{array}{c}6\\ 5\\ 4\\ 3\\ 2\end{array}\right]\Vert }_{\infty }={\Vert \left[\begin{array}{c}12\\ 15\\ 13\\ 11\\ 9\end{array}\right]\Vert }_{\infty }=15,$$

Full text

I Introduction

Automated controller synthesis for systems subject to an a priori given specification is an attractive means of controller design; such correct-by-construction synthesis methods have attracted considerable interest in the past years [1]. However, all methods that are capable of solving general problems face fundamental limitations in terms of the system dimensionality and specification complexity that can be handled. A large body of work resorts to abstraction procedures [2], where a finite representation of the system is obtained for which synthesis techniques for finite-state systems can be applied, but naive grid-based abstraction methods quickly run into the curse of dimensionality since, in general, the number of discrete states grows exponentially with system dimension. Other methods operate directly on a continuous state space but also “blow up” when the representation and treatment of high-dimensional sets become overwhelmingly expensive.

Ways to overcome these limitations in the abstraction step have been proposed. A large problem can sometimes be decomposed into several lower-dimensional problems [3, 4], where each subproblem is easier to solve. Such methods typically require the decomposition to be given, and that the coupling between the subproblems is relatively weak. Similarly, compositional abstraction procedures have been suggested as a method to construct abstractions with a size that does not scale exponentially with dimensionality [5, 6]; to capture the dynamics more accurately the abstract subsystems can be made partly overlapping [7]. These compositional approaches account for dynamic coupling between subsystems; as a contrast, we consider in this paper subsystems that are dynamically independent but coupled through a common specification. Other work has exploited monotonicity [8, 9] to alleviate the curse of dimensionality when constructing an abstraction, and it has been found that under certain stability conditions discretization can be circumvented completely by constructing an abstraction whose states are a finite-horizon input history [10, 11].

In this work we explore a different approach to tackle high-dimensionality: exploitation of symmetries. The system we consider is the aggregate system consisting of $N$ almost-homogeneous switched subsystems, where $N$ can be very large. Due to homogeneity the overall system exhibits symmetries in the dynamics which can be leveraged to enable synthesis for a special type of symmetric constraints called counting constraints. Like many existing methods in the literature, our solution approach is abstraction-based, but with the crucial difference that we construct only one abstraction that is used for all subsystems, thus roughly reducing the number of states in the abstraction from $\mathcal{O}(1/\eta^{Nn_{x}})$ in the naïve approach to $\mathcal{O}(1/\eta^{n_{x}})$, where $n_{x}$ is the dimension of a single subsystem and $\eta$ the precision. Our main contributions are not abstraction construction itself, but showing how a single abstraction can be utilized for control synthesis purposes.111We consider here the problem of synthesizing open-loop infinite-horizon trajectories that satisfy counting constraints for a given initial condition, rather than the more general problem of synthesizing feedback controllers that enforce specifications for sets of initial conditions. We are able to account for mild heterogeneity among the subsystems by assuming certain stability conditions. Utilizing a single abstraction for multiple subsystems allows us to synthesize controlled trajectories for problems of very high dimension that exhibit these symmetries—we demonstrate the method on a $10,000$-dimensional switched system with $2^{10,000}$ modes. If there are several distinct classes of subsystems, one abstraction can be constructed for each class; and complexity increases linearly with the number of classes.

This work is motivated by the problem of scheduling of thermostatically controlled loads. Examples of TCLs include air conditioners, water heaters, refrigerators, etc., that operate around a temperature set point by switching between being on and off. TCL owners are typically indifferent to small temperature perturbations and accept temperatures in a range around their desired set point; this range is called the dead band. The idea behind TCL scheduling is that an electric utility company can leverage the implied flexibility—–which becomes meaningful for large collections of TCLs–—to shape aggregate demand on the grid. Previous work has resulted in control algorithms based on broadcasting a universal set-point temperature [12], Markov chain bin models [13, 14], and priority stacks [15], with the objective to track a power signal. However, they do not take into account any hard constraints on aggregate power consumption (i.e., the number of TCLs in mode on), which is crucial to do to avoid overloading the grid or underutilizing the power generated by renewable resources.

The fundamental problem in TCL scheduling is to simultaneously meet local safety constraints (i.e., maintain each TCL in its deadband), and global aggregate constraints. These constraints are special instances of counting constraints, which is the type of constraints we consider in this work. The TCL scheduling example will be used throughout the paper to motivate our discussion:

Example 1**.**

Let $\{\mathbf{x}_{n}\}_{n\in[N]}$ be the states of a family of $N$ switched systems of the following form:

[TABLE]

where $f_{\texttt{on}}$ and $f_{\texttt{off}}$ are functions $\mathbb{R}\times\mathcal{D}\rightarrow\mathbb{R}$ for some bounded disturbance set $\mathcal{D}$. For a desired dead band $[\underline{a},\overline{a}]$, and bounds $[\underline{K},\overline{K}]$ on the aggregate number of TCLs that are in mode on, synthesize switching protocols that enforce

[TABLE]

Our approach can be seen as a marriage between “bin” abstraction previously used in the TCL literature, and schedule search via integer programming previously considered e.g. for intersection controllers [16]—generalized to arbitrary incrementally stable systems and an infinite horizon.

A preliminary version of this work appeared in [17] where the mode-counting problem was introduced for homogeneous collections of systems. In this paper we generalize the counting problem to encompass both mode- and state-counting constraints; make the construction of the abstractions robust to model deviations—thus allowing mild heterogeneity; improve the rounding scheme that enables search for solutions via (non-integer) linear programming; and sharpen analytical results. A different extension of counting problems that allows for a richer class of specifications—temporal logic formulas defined over counting propositions—was reported in [18].

I-A Overview and Paper Structure

We give a high-level overview of the proposed solution approach to the problem of synthesizing controllers for large collections of systems with counting constraints, namely the counting problem as formally defined in Section III. By exploiting the fact that the collection is almost homogeneous, we first construct a single finite abstraction that approximates all the systems in the collection (Section IV). By means of this abstraction, we transform a continuous-state counting problem into a discrete one. We then treat this finite abstraction as a histogram where each abstract state corresponds to a bin of the histogram and counts the number of systems whose continuous-state is in some equivalence class associated with that abstract state (Section V). As the continuous states evolve in time according to the chosen control inputs, the counts on this histogram also evolve according to certain aggregate dynamics. In particular, the evolution of the histogram can be represented as a constrained linear system over an integer lattice, and the counting constraints become linear constraints on the states and input of the aggregate system (Section V-A). This representation allows us to reduce the discrete counting problem to the feasibility of an integer linear program (Section V-C). Through this integer linear program, we search for inputs to steer the systems to periodic trajectories, captured by the cycles of the abstraction graph, so that the satisfaction of the counting constraints can be guaranteed indefinitely. We then analyze this integer linear program and a relaxation of the integrality constraints, and their relation to the (in)-feasibility of the original problem (Section VI). An extension of the theory to the multi-class setting is briefly presented in Section VII. To highlight the applicability of the method, one numerical example and one application-motivated example are provided in Section VIII before the paper is summarized in Section IX. In an effort to improve readability, technical proofs are deferred to the appendix.

II Notation and Preliminaries

We introduce some notation that is used throughout the paper. The set of real numbers is denoted $\mathbb{R}$, the set of positive reals $\mathbb{R}_{+}$, and the set of non-negative integers $\mathbb{N}$. To express a finite set of positive integers, we write $[N]=\{0,\ldots,N-1\}$. The indicator function of a set $X$ is denoted $\mathds{1}_{X}(x)$ and is equal to 1 if $x\in X$ and to [math] otherwise. The identity function on a set $A$ is written as $\text{Id}_{A}$. For two sets $X$ and $Y$, we write the Minkowski sum as $X\oplus Y=\{x+y:x\in X,y\in Y\}$, and the Minkowski difference $X\ominus Y=\{x:\{x\}\oplus Y\subset X\}$. Set complement is denoted $X^{C}$.

To denote the floor and ceiling of a number we write $\left\lfloor\cdot\right\rfloor$ and $\left\lceil\cdot\right\rceil$. We use the same notation for vectors, where the operations are performed component-wise. We write the infinity norm as $\|\cdot\|_{\infty}$ and the $1$-norm as $\|\cdot\|_{1}$. The $\epsilon$-ball in $p-$norm centered at the point $x$ is denoted as $\mathcal{B}_{p}(x,\epsilon)=\{y:\|y-x\|_{p}\leq\epsilon\}$. The vector of all $1$’s is written as $\mathbf{1}$. For a function $f:X\rightarrow\mathbb{R}$ with a finite domain we abuse notation and write $\|f\|_{1}=\sum_{x\in X}|f(x)|$ for the “1-norm” of the finite image set. The function that is constantly equal to [math] is written $\mathbf{0}$.

Given an ODE $\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{x}(t)=f_{\mu}(\mathbf{x}(t),\mathbf{d}(t))$, where $\mathbf{d}(t)$ is an uncontrolled input, the corresponding flow operator is denoted $\phi_{\mu}(t,x,\mathbf{d})$ and has the properties that $\phi_{\mu}(0,x,\mathbf{d})=x$, $\frac{\mathrm{d}}{\mathrm{d}t}\phi_{\mu}(t,x,\mathbf{d})=f_{\mu}(\phi_{\mu}(t,x,\mathbf{d}),\mathbf{d}(t))$. A $\mathcal{KL}$-function $\beta:\mathbb{R}_{+}\times\mathbb{R}_{+}\rightarrow\mathbb{R}_{+}$ is characterized by being strictly increasing from 0 in the first argument and decreasingly converging to 0 in the second argument.

II-A Transition Systems and Simulations

We employ the following definition for a transition system, which captures systems with both continuous and discrete state spaces.

Definition 1**.**

A transition system (TS) is a tuple $\Sigma=(\mathcal{Q},\mathcal{U},\longrightarrow,Y)$, where $\mathcal{Q}$ is a set of states, $\mathcal{U}$ a set of actions (inputs), $\longrightarrow\subset\mathcal{Q}\times\mathcal{U}\times\mathcal{Q}$ a transition relation, and $Y:\mathcal{Q}\longrightarrow\mathcal{Y}$ an output function. We say that $\Sigma$ is a deterministic finite transition system (DFTS) if i) transitions are deterministic, i.e., if $(q,\mu,q^{\prime})\in\mathcal{\longrightarrow}$ and $(q,\mu,q^{\prime\prime})\in\mathcal{\longrightarrow}$, then $q^{\prime}=q^{\prime\prime}$, and ii) $\mathcal{Q}$ is finite.

For simplicity, existence of a transition $(q,\mu,q^{\prime})\in\longrightarrow$ is written $q\overset{\mu}{\longrightarrow}q^{\prime}$. By a trajectory of a transition system, we mean a sequence $\mathbf{x}(0)\mathbf{x}(1)\mathbf{x}(2)\ldots$ of states in $\mathcal{Q}$ with the property that $\mathbf{x}(s)\overset{\sigma(s)}{\longrightarrow}\mathbf{x}(s+1)$ for some $\sigma(s)\in\mathcal{U}$, for all $s\in\mathbb{N}$.

In this paper we restrict attention to finite input sets $\mathcal{U}$; such transition systems are obtained when a continuous-time switched system is time-discretized, and in this case the inputs are the modes of the switched system. A switching protocol is then a function $\pi$ that generates control inputs from information about the current state222A switching protocol may also have internal memory states.. A trajectory generated by a switching protocol $\pi$ is a trajectory where $\sigma(s)$ is generated by $\pi$.

For two systems with the same input space $\mathcal{U}$ and normed output space $\mathcal{Y}$, we adopt the following notion of system bisimilarity [19].

Definition 2**.**

Two transition systems $(\mathcal{Q}_{1},\mathcal{U},\longrightarrow_{1},Y_{1})$ and $(\mathcal{Q}_{2},\mathcal{U},\longrightarrow_{2},Y_{2})$ are $\epsilon$-approximately bisimilar if there exists a relation $R\subset\mathcal{Q}_{1}\times\mathcal{Q}_{2}$ such that the sets $R(q_{1})=\{q_{2}:(q_{1},q_{2})\in R\}$ and $R^{-1}(q_{2})=\{q_{1}:(q_{1},q_{2})\in R\}$ are non-empty for all $q_{1},q_{2}$, and such that for all $(q_{1},q_{2})\in R$,

$\|Y_{1}(q_{1})-Y_{2}(q_{2})\|_{\infty}\leq\epsilon$, 2. 2.

if $q_{1}\overset{\mu}{\longrightarrow_{1}}q_{1}^{\prime}$, there exists $q_{2}\overset{\mu}{\longrightarrow_{2}}q_{2}^{\prime}$ s.t $(q_{1}^{\prime},q_{2}^{\prime})\in R$, 3. 3.

if $q_{2}\overset{\mu}{\longrightarrow_{2}}q_{2}^{\prime}$, there exists $q_{1}\overset{\mu}{\longrightarrow_{1}}q_{1}^{\prime}$ s.t $(q_{1}^{\prime},q_{2}^{\prime})\in R$.

II-B Graphs

The following standard notions are used for a directed graph $G=(\mathcal{Q},E)$ with node set $\mathcal{Q}$ and edge set $E\subset\mathcal{Q}\times\mathcal{Q}$. A path in $G$ is a list of edges $(q_{0},q_{1})(q_{1},q_{2})\ldots(q_{J-1},q_{J})$. The distance between two nodes $q_{0}$ and $q_{J}$ is the length of the shortest path connecting the two. The graph diameter $\operatorname{diam}(G)$ is the longest distance between two nodes in the graph. If the first and last nodes in a path are equal, i.e. $q_{J}=q_{0}$, the path is a cycle. A cycle is simple if it visits every node at most one time. For a subset of nodes $D\subset\mathcal{Q}$ (or the corresponding subgraph, which we use interchangeably), it is said to be strongly connected if for each node pair $(q,q^{\prime})\in D$, there exists a path from $q$ to $q^{\prime}$. Any directed graph can be decomposed into strongly connected components. The period of a subgraph $D$ is the greatest common divisor of the lengths of all cycles in $D$. A subgraph $D$ is called aperiodic if it has period one.

III The Counting Problem

We first define the concept of a counting constraint.

Definition 3**.**

A counting constraint $(X,R)$ for a collection of $N$ subsystems numbered from [math] to $N-1$, with the same state space $\mathcal{Q}$ and input space $\mathcal{U}$, is a set $X=X^{\mathcal{Q}}\times X^{\mathcal{U}}$ with $X^{\mathcal{Q}}\subset\mathcal{Q}$ and $X^{\mathcal{U}}\subset\mathcal{U}$, and a bound $R$. We say that the counting constraint is satisfied by state-input pairs $(q_{n},\mu_{n})$, $n\in[N]$, if the number of pairs that fall in $X$ is less than or equal to $R$:

[TABLE]

Remark 1**.**

For simplicity of exposition we restrict attention to counting sets in the form of a product of a set $X^{\mathcal{Q}}$ in state space and a set $X^{\mathcal{U}}$ in mode space. However, our solution method in Section V can handle more general subsets of $\mathcal{Q}\times\mathcal{U}$.

This notion of constraint extends that in [17] by allowing $X$ to be a subset of $\mathcal{Q}\times\mathcal{U}$, instead of restricting it to a singleton subset of $\mathcal{U}$. This generalization permits counting constraints that include the state space, as opposed to only counting the number of systems using each input (i.e., mode-counting).

Example 2**.**

The TCL scheduling constraints (2) and (3) fit into this class. Let the output space be equal to the state space, i.e. $\mathcal{Y}=\mathbb{R}$ and $Y(x)=x$. Then the following counting constraints specify that the number of TCLs that are in mode on at time $t$ should be in the interval $[\underline{K},\overline{K}]$:

[TABLE]

In addition, the dead band constraints (2) can be imposed by the counting constraint

[TABLE]

The goal of this work is to synthesize controllers that satisfy counting constraints for a collection of switched systems with states $\mathbf{x}_{i}(t)$, $i\in[N]$, governed by dynamics

[TABLE]

for $\mathbf{x}_{i}(t)\in\mathbb{R}^{n_{x}}$ and a disturbance signal $\mathbf{d}_{i}(t)$ assumed to take values in a set $\mathcal{D}$. Although the system description is identical for all $N$ systems in the family, the disturbance can model mild heterogeneity such as modeling inaccuracies and parameter variations (e.g., $\mathbf{d}_{i}(t)=\mathbf{d}_{i})$ across subsystems.

We assume that the vector fields and the disturbance signals satisfy standard assumptions for existence and uniqueness of solutions. In addition, we assume that the vector fields exhibit a certain form of stability [20].

Assumption 1**.**

For each $\mu\in\mathcal{U}$, the vector field $f_{\mu}(x,d)$ is $C^{1}$ in $x$ and continuous in $d$. Furthermore, the nominal system $\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{x}=f_{\mu}(\mathbf{x},\mathbf{0})$ is forward-complete and incrementally stable. That is, there exists a $\mathcal{KL}$-function $\beta_{\mu}$ such that

[TABLE]

Note that the vector field $f_{\mu}(x,d)$ being $C^{1}$ in $x$ implies that $f_{\mu}(x,d)$ is Lipschitz in $x$ (with some Lipschitz constant $K_{\mu}$) when the states are constrained to a compact set.

Furthermore, we assume that the disturbance signal is continuous, and that its effect is bounded in an absolute sense.

Assumption 2**.**

The disturbance signal $\mathbf{d}:\mathbb{R}\rightarrow\mathcal{D}$ is continuous. Furthermore for all $\mu\in\mathcal{U}$, compared to the nominal vector fields without disturbance, the effect of the disturbance is less than $\bar{\delta}_{\mu}$:

[TABLE]

To leverage the notion of bisimilarity we consider the time-sampled counterpart of (6), and define the problem we want to solve in terms of the resulting transition system rather than the ODE. An ODE solution where state counting constraints are violated in between samplings (but satisfied at sample instants) is still a valid solution to a time-sampled counting problem. If such inter-sample violations are unacceptable, the counting sets can be contracted by some margin determined by the dynamics to ensure satisfaction for all $t\in\mathbb{R}_{+}$ [21].

We proceed with a problem definition for transition systems. Consider a family of $N$ identical subsystems with dynamics given by a transition system and some initial conditions $\{\mathbf{x}_{n}(0)\}_{n\in[N]}$. The problem we seek to solve is the following:

Problem 1**.**

Consider $N$ subsystems, all governed by the same transition system $\Sigma_{TS}=(\mathcal{Q},\mathcal{U},\longrightarrow,Y)$, where $\mathcal{U}$ is finite, and assume that initial states $\{\mathbf{x}_{n}(0)\}_{n\in[N]}$ and $L$ counting constraints $\{X_{l},R_{l}\}_{l\in[L]}$ are given. Synthesize individual switching protocols $\{\pi_{n}\}_{n\in[N]}$ such that the generated actions $\sigma_{n}(0)\sigma_{n}(1)\sigma_{n}(2)\ldots$ and trajectories $\mathbf{x}_{n}(0)\mathbf{x}_{n}(1)\mathbf{x}_{n}(2)\ldots$ for $n\in[N]$ satisfy the counting constraints

[TABLE]

Let $(N,\Sigma_{TS},\{\mathbf{x}_{n}(0)\}_{n\in[N]},\{X_{l},R_{l}\}_{l\in[L]})$ denote an instance of this problem.

Naïvely, if $\mathcal{Q}$ is a $n_{q}$-dimensional space, the overall system can be viewed as a $N\times n_{q}$-dimensional hybrid system with $|\mathcal{U}|^{N}$ modes and is thus beyond the reach of standard synthesis techniques for large $N$. However, from a specification point of view it does not matter which subsystems that contribute to the summation in (8), so Problem 1 exhibits subsystem permutation symmetry in its specification. In addition, there is also subsystem permutation symmetry in the dynamics since all subsystem states are governed by the same transition system. As an implication there is no need to keep track of identities of individual subsystems—they are all equivalent from both a dynamics- as well as from a specification-point of view. This symmetry (i.e., permutation invariance) is the key to solving large-scale counting problems.

In the following, we construct the time-sampled analogue of (6) in order to leverage the notion of bisimilarity. We also restrict the domain to a bounded set $\mathcal{X}\subset\mathbb{R}^{n_{x}}$ and let the $\tau$-sampled counterpart of (6) confined to $\mathcal{X}$ be

[TABLE]

where $x\overset{\mu}{\underset{\tau}{\longrightarrow}}x^{\prime}$ if and only if there exists $\mathbf{d}:[0,\tau]\rightarrow\mathcal{D}$ such that $x^{\prime}=\phi_{\mu}(\tau,x,\mathbf{d})$. The instance

[TABLE]

of Problem 1 is referred to as the continuous-state counting problem. Since $\Sigma_{\tau}$ is a non-deterministic, infinite transition system, this is a difficult problem to solve. In the following we construct a DFTS $\Sigma_{\tau,\eta}$ that is approximately bisimilar to $\Sigma_{\tau}$, and relate the corresponding problem instances.

IV Abstracting the Continuous-State Counting Problem

We first recall the abstraction procedure from [22] to create a finite-state model of (6). Under the assumptions outlined above, we establish approximate bisimilarity between the continuous-state system and its finite-state abstraction, which enables us to state results in Section IV-B relating the solvability of the corresponding counting problems.

IV-A Abstraction Procedure

Consider a system of the form (6). For a state discretization parameter $\eta>0$ we define an abstraction function $\kappa_{\eta}:\mathcal{X}\mapsto\mathcal{X}$ as $\kappa_{\eta}(x)=\eta\cdot\left\lfloor x/\eta\right\rfloor+(\eta/2)\mathbf{1}$. This function is constant on hyperboxes of side $\eta$, and its image of the compact set $\mathcal{X}$ is finite. The abstraction function defines the $(\tau,\eta)$-discretized counterpart of (6) as the transition system

[TABLE]

where $q\overset{\mu}{\underset{\tau,\eta}{\longrightarrow}}q^{\prime}$ if and only if $\kappa_{\eta}(\phi_{\mu}(\tau,q,\mathbf{0}))=q^{\prime}$.

The domain $\mathcal{X}$ is partitioned into uniform boxes of size $\eta$ that represent discrete states, and transitions are established by simulating each mode, without disturbance, during a time $\tau$, starting at the center of the boxes. As opposed to $\Sigma_{\tau}$, the resulting transition system $\Sigma_{\tau,\eta}$ is finite and deterministic—for each state $q$ and action $\mu$ there exists (at most) one successor state $q^{\prime}$. This makes the resulting abstract counting problem easier to solve. Similarly to a result from [22] we can now show that an abstraction constructed in this way is bisimilar to the time-sampled system (9) if a certain inequality holds.

Proposition 1**.**

Assume that Assumptions 1 and 2 hold, i.e., for all $\mu\in\mathcal{U}$ there are Lipschitz constants $K_{\mu}$, $\mathcal{KL}$-functions $\beta_{\mu}$, and disturbance effect bounds $\overline{\delta}_{\mu}$ associated with the modes of (6). Then, if for all $\mu\in\mathcal{U}$,

[TABLE]

the $(\tau,\eta)$-discretized abstraction $\Sigma_{\tau,\eta}$ and the $\tau$-sampled system $\Sigma_{\tau}$ are $\epsilon$-approximately bisimilar.

The trajectories of $\epsilon$-approximately bisimilar systems remain within distance $\epsilon$ of each other [19], and the additional term $\frac{\overline{\delta}_{\mu}}{K_{\mu}}\left(e^{K_{\mu}\tau}-1\right)$ in (11) introduces a robustness margin that allows the same abstraction to be bisimilar to a family of mildly heterogeneous subsystems. Using these facts we can establish relations between existence of solutions of the counting problem in the continuous-state and discrete-state settings.

IV-B Relations Between the Continuous-State and Discrete-State Counting Problems

If $\epsilon$-approximate bisimilarity holds, trajectories of $\Sigma_{\tau,\eta}$ and $\Sigma_{\tau}$ are guaranteed to remain $\epsilon$-close when initial states are $\epsilon$-close and corresponding actions are chosen. Therefore, solvability of the counting problem is equivalent for the two up to an approximation margin $\epsilon$. To make precise statements we introduce functions $\mathcal{G}^{\pm\epsilon}$ that expand (resp. contract) a counting set $X=X^{\mathcal{X}}\times X^{\mathcal{U}}$ in $\mathcal{X}$-space before quantization:

[TABLE]

Now we can describe how solutions of the counting problem can be mapped between the time-sampled system $\Sigma_{\tau}$ and the time-sampled and state-quantized system $\Sigma_{\tau,\eta}$.

Theorem 1**.**

Let $\Sigma_{\tau}$ and $\Sigma_{\tau,\eta}$ be the time-sampled and time-sampled and state-quantized systems constructed from a system on the form (6), such that $\Sigma_{\tau}$ and $\Sigma_{\tau,\eta}$ are $\epsilon$-approximately bisimilar. Let $\kappa_{\eta}$ be the abstraction function for $\Sigma_{\tau,\eta}$.

If there exists a solution to the instance

[TABLE]

of Problem 1, then there exists a solution to the instance

[TABLE]

An converse result can also be obtained, with the only difference that an additional correction $\eta/2$ is required to account for counting sets that are inflated by quantization.

Theorem 2**.**

Under the same assumptions as in Theorem 1, if there is no solution to the instance

[TABLE]

of Problem 1, then there is no solution to the instance

[TABLE]

Example 3** (continued).**

We proceed with the TCL scheduling problem by adjusting the constraints. The mode-counting constraints (4) lack a state space part, and are therefore not modified. However, by Theorem 1 the local safety constraints (5) must be expanded to

[TABLE]

in order for a solution of the discrete-state instance to be mappable to a valid solution of the continuous-state instance of Problem 1. Similarly, by Theorem 2, if it can be shown that the discrete-state instance obtained by contracting the counting set to $[\underline{a}-\epsilon-\eta/2,\overline{a}+\epsilon+\eta/2]^{C}$ before discretizing is infeasible, then also the original problem is infeasible.

V Solving the Discrete-State Counting Problem

Having described the reduction of a continuous-state instance of Problem 1 to a deterministic discrete-state instance, we proceed with a solution procedure for the latter. Consider a DFTS $\Sigma_{DFTS}=(\mathcal{Q},\mathcal{U},\longrightarrow,Y)$; it may be the result of a continuous-state abstraction or just represent a discrete structure onto which a counting problem is defined. A DFTS can alternatively be viewed as a directed graph $G=(\mathcal{Q},\longrightarrow)$, where $\mathcal{Q}$ is a set of nodes and $\longrightarrow$ is a set of edges, and each edge is labeled with an action from $\mathcal{U}$. This dual viewpoint allows us to leverage graph-theoretical concepts to investigate the aggregate system; in the following we use the transition system viewpoint and the graph viewpoint interchangeably. We first define aggregate dynamics which take the form of a constrained linear system. The solution method is then presented in the form of a linear feasibility problem over aggregate states.

V-A Aggregate Dynamics as a Linear System

Consider a total of $N$ subsystems whose dynamics are governed by a DFTS $\Sigma_{DFTS}=(\mathcal{Q},\mathcal{U},\longrightarrow,Y)$.

We introduce $|\mathcal{Q}|$ aggregate states labeled $w_{q}$ for $q\in\mathcal{Q}$, that describe the number of individual systems that are at state $q$. By also introducing control actions $r_{q}^{\mu}$ that represent the number of systems at state $q$ using action $\mu$, the aggregate dynamics can be written as

[TABLE]

where $\mathcal{N}_{q}^{\mu}=\left\{q^{\prime}\in\mathcal{Q}:q^{\prime}\overset{\mu}{\longrightarrow}q\right\}$ is the set of predecessors of $q$ under the action $\mu$. We constrain the control actions $r_{q}^{\mu}$ such that for all $\mu\in\mathcal{U}$,

[TABLE]

which ensures the continued positivity of the states: $w_{q}(s+1)\geq 0$ for $q\in\mathcal{Q}$. Furthermore, the invariant $\|\mathbf{w}(s)\|_{1}=N$ holds over time, where $N$ is the total number of subsystems. In the following, we use the compact notation

[TABLE]

to denote this system, where $B$ is composed of the incidence matrices $A^{\mu},\mu\in\mathcal{U}$ of the system graph. The state space $\mathcal{W}$ and (state-dependent) admissible control set $\mathcal{R}$ of this system are then

[TABLE]

When the discrete-state counting problem is obtained via the abstraction procedure in Section IV, the aggregate dynamics (16) can be seen as a time-varying histogram, where each bin in the histogram represents the number of subsystems in a given area of the state space, as illustrated in Fig. 1.

Counting constraints become linear constraints on these aggregate states, so our problem can be rephrased as a linear feasibility problem. However, the feasibility problem is infinite-dimensional since we are searching for solution-trajectories over an infinite horizon. To limit the search to a finite-dimensional space we consider solutions that consist of a finite prefix trajectory and a periodic suffix trajectory. While the prefix trajectory is just a trajectory of the aggregate dynamics, we consider suffix trajectories that are defined in terms of cycle assignments.

V-B Cycle assignments

We consider cycles $C$ of the form

[TABLE]

of length given by $|C|=I$, where $(q_{i},\mu_{i},q_{{i+1}})\in\longrightarrow$. A cycle assignment assigns “weights” (or numbers) of subsystems to each node along a cycle. If the graph represents an abstraction, the assignment counts the number of subsystems whose continuous states are in the vicinity of the abstract states on this cycle.

Definition 4**.**

An assignment to a cycle $C$ is a function $\alpha:[|C|]\rightarrow\mathbb{R}_{+}$. If $\alpha(i)\in\mathbb{N}$ for all $i\in[|C|]$, $\alpha$ is an integer assignment.

For now we will only be concerned with integer assignments that represent a number of subsystems. However, later in the paper we will study relaxed non-integer solutions consisting of assignments that are not necessarily integral.

An assignment assigns subsystems to a cycle, so that the state of a subsystem circulates along the abstract states on this cycle as time progresses (provided that the appropriate actions are chosen). The movement corresponds to a circular shift of the assignment.

Definition 5**.**

For an assignment $\alpha$ we define its $s$-step circulation, denoted $\alpha^{\circlearrowleft s}:[|C|]\rightarrow\mathbb{R}_{+}$, as the shifted function

[TABLE]

The periodicity is manifested by the relation $\alpha^{\circlearrowleft(|C|+s)}=\alpha^{\circlearrowleft s}$. To capture how counting quantities vary during the circulation we introduce the following notation for the matching of a cycle and a circulated assignment.

Definition 6**.**

For a cycle $C$ and an assignment $\alpha:[|C|]\rightarrow\mathbb{R}_{+}$ the $X$-count at time $s$ of a counting set $X$ is defined as

[TABLE]

If the cycle $C$ includes elements from a set $X$ coming from a counting constraint, (20) counts the number of subsystems contributing to the counting constraint at time $s$, assuming that the subsystems make up the assignment $\alpha$ at time zero and circulate in the cycle $C$. The concept is illustrated in Fig. 2.

Next we introduce a function that returns the highest count over all circulations for given cycle-assignment pairs.

Definition 7**.**

The maximal $X$-count for a cycle $C$ with assignment $\alpha$, denoted $\operatorname{maxcnt}^{X}(C,\alpha)$, is the maximal number of subsystems simultaneously in $X$ when $\alpha$ is circulated around $C$:

[TABLE]

The maximal $X$-count can be computed as the maximum entry in a matrix-vector product $\operatorname{maxcnt}^{X}(C,\alpha)=\left\|B_{C}^{X}\alpha\right\|_{\infty}$, where $B_{C}^{X}\in\{0,1\}^{|C|\times|C|}$ is a circulant binary matrix.

Example 4**.**

Going back to cycle and assignment in Fig. 2, it can be seen that the maximal $X$-count occurs for $s=1$, and that $\operatorname{maxcnt}^{X}(C,\alpha)=15$. The same value is obtained as the infinity norm of the matrix-vector product

[TABLE]

where the first row of the matrix has ones at positions in $C$ corresponding to states in $X$ (i.e. $q_{1}$, $q_{2}$, and $q_{3}$), and the remaining rows are shifted to form a circulant matrix.

Our solution method to the counting problems requires consideration of multiple assignments that circulate in different cycles. For this reason we introduce the joint maximal count of assignments to cycles.

Definition 8**.**

The maximal joint $X$-count for a set of cycles $\{C_{j}\}_{j\in J}$ and a matching set of assignments $\{\alpha_{j}\}_{j\in J}$, denoted $\operatorname{maxcnt}^{X}\left(\{C_{j}\}_{j\in J},\{\alpha_{j}\}_{j\in J}\right)$, is the maximal number of subsystems simultaneously in $X$ when the assignments $\{\alpha_{j}\}_{j\in J}$ are synchronously circulated around the cycles $\{C_{j}\}_{j\in J}$:

[TABLE]

Also the maximal joint $X$-count can be expressed as the maximal element in a sum of matrix-vector products.

Proposition 2**.**

The maximal joint $X$-count satisfies

[TABLE]

where $k_{j}=\operatorname{lcm}\left(\{|C_{j}|\}_{j\in J}\right)/|C_{j}|$, $\operatorname{lcm}$ is the least common multiple, $\mathbf{1}_{k}$ is the length $k$ column vector of all ones, and $\otimes$ is the Kronecker product.

The matrices in (21) are large if the least common multiple is large, which illustrates that relative assignment circulations must be taken into account. We elaborate further on these matters, and how they are connected to graph periodicity, in Section V-D below.

V-C Solution to the Discrete-State Counting Problem

Consider now a problem instance

[TABLE]

and let $\mathbf{w}(0)$ be the aggregate initial state

[TABLE]

We restrict attention to solutions of a particular form; it is shown in later sections that this is without loss of generality.

Definition 9**.**

A trajectory for an aggregate initial state $\mathbf{w}(0)$ is of prefix-suffix type if it consists of a finite number of inputs $\mathbf{r}(0),\ldots,\mathbf{r}(T-1)$, and a set of cycles $\{C_{j}\}_{j\in J}$ with assignments $\{\alpha_{j}\}_{j\in J}$ such that the cycles are populated with their respective cycle assignments at time $T$.

Given aggregate initial states $\mathbf{w}(0)$, a set $\{C_{j}\}_{j\in[J]}$ of cycles in $G$, and a prefix horizon $T$, prefix-suffix trajectories can be extracted from feasible points of the following linear feasibility problem:

[TABLE]

The constraints (24a) and (24b) enforce counting constraints in the prefix and suffix phases, and (24c) ensures that the prefix and suffix phases are adequately connected. The remaining constraints (24d)-(24f) certify that the prefix trajectory is dynamically feasible with respect to (16).

The number of variables and (in)equalities in (24) are $\mathcal{O}(T|\mathcal{Q}||\mathcal{U}|+\sum_{j\in[J]}|C_{j}|)$ and $\mathcal{O}\left(LT+L\operatorname{lcm}\left(\{|C_{j}|\}_{j\in J}\right)+T|\mathcal{Q}|\right)$, respectively (not counting positivity constraints that solvers handle easily). Crucially, these numbers do not depend on the total number of subsystems $N$ which makes the approach suitable for large $N$ and moderate-sized graphs.

The set $\{C_{j}\}_{j\in J}$ of cycles is an input to the optimization problem (24); the solver selects good cycles within this set by setting the assignments of the rest to zero. Section VI presents results regarding when an input cycle set is “sufficiently rich” not to compromise the existence of solutions. On the other hand, in practice, we find that randomly sampling a set of cycles is sufficient for feasibility in many problems.

V-D Size Reductions of the Linear Program

As we point out next, there are certain ways to further decrease the number of variables and/or inequalities—sometimes without loss of generality.

First of all, state-mode pairs in $G$ with mandated [math]-count can be pruned from $G$ to decrease the number of aggregate variables. Specifically, we can construct a pruned graph $\tilde{G}$ in the following way: if $(q,\mu)\in X_{l}$ and $R_{l}=0$, then remove the action $\mu$ at $q$. If this results in nodes in $\tilde{G}$ with no outgoing edges, prune these nodes together with incoming edges and repeat until all nodes have at least one valid action. This procedure is equivalent to finding the largest controlled-invariant set contained in $X_{l}^{C}$.

Secondly, the following result drastically reduces the number of suffix counting constraint inequalities in (24b).

Proposition 3**.**

The joint $X$-count for two cycles $C_{0}$ and $C_{1}$ with co-prime length, i.e. $\gcd(|C_{0}|,|C_{1}|)=1$, can be computed as

[TABLE]

Proof.

If $|C_{0}|$ and $|C_{1}|$ are coprime, by the Chinese Remainder Theorem [23], the equations $k_{1}=s\mod|C_{1}|$ and $k_{2}=s\mod|C_{2}|$ have a unique solution $s<|C_{0}||C_{1}|$ for every pair $k_{1}<|C_{1}|,k_{2}<|C_{2}|$. It follows that every circulation of $\alpha_{0}$ in $C_{0}$ at some point coincides with every circulation of $\alpha_{1}$ in $C_{1}$; hence the maximal joint $X$-count is equal to the sum of individual maximal $X$-counts for the two cycles. ∎

The proof illustrates that if there is no periodicity, every relative assignment position will eventually be attained, including the worst-case relative assignment position which is exactly the combined worst-case absolute positions of individual cycles. However, in the presence of periodicity only a subset of all relative positions will be attained and representing the worst-case within this subset is not as straight-forward.

More generally, if the cycle set can be partitioned into sets of cycles with mutually co-prime length, the number of inequalities can be reduced. That is, for sets of cycles $\left\{C_{j}\right\}_{j\in J}$ and $\left\{C_{j^{\prime}}\right\}_{j^{\prime}\in J^{\prime}}$ with the property that $\text{gcd}(|C_{j}|,|C_{j^{\prime}}|)=1$ for all pairs $j\in J$, $j^{\prime}\in J^{\prime}$, it holds that

[TABLE]

Example 5**.**

To exemplify this reduction, consider a set of cycles with lengths ranging from 2 to 20. We have $\operatorname{lcm}([21])=232792560$ and $\operatorname{lcm}([21]\setminus\{11,13,17,19\})+11+13+17+19=5100$. The number 232792560 is the number of inequalities in the naïve approach (21), which by (25) can be drastically reduced to 5100 if the cycles of prime lengths $11,13,17$, and $19$ are considered separately.

If the number of constraints is still prohibitively large, it can be replaced by a conservative constraint as follows.

Remark 2**.**

The constraint (24b) can be substituted by the conservative constraint

[TABLE]

which groups cycles by cycle length $i$ and disregards effects from periodicity. The number of constraints in this case becomes $L(1+\|\bigcup_{j\in J}\{|C_{j}|\})\|_{1})$ instead of $L\operatorname{lcm}\left(\{|C_{j}|\}_{j\in J}\right)$, where $\bigcup_{j\in J}\{|C_{j}|\}$ is the set of cycle lengths without repetition.

While the total number of subsystems $N$ does not impact the number of inequalities or constraints, it might affect the performance of integer linear program solvers since the number of candidate integer points grows with $N$. In addition, the converse results in Section VI depend on $N$—in order to prove infeasiblity of the problem very large horizons $T$ and/or cycle sets may be required. If $N$ is prohibitively large for this purpose it can be scaled down if a certain divisibility condition holds: if there is a common divisor $S$ that divides $w_{q}(0)$ for all $q\in\mathcal{Q}$, and that divides $R_{l}$ for all $l\in[L]$, then there is a 1-1 correspondence between solutions of (24) and solutions of its analogue obtained from the substitutions $\mathbf{w}(0)\rightarrow\mathbf{w}(0)/S$ and $R_{l}\rightarrow R_{l}/S$. The correspondence simply consists in scaling $\mathbf{r}$, $\mathbf{w}$, and the assignments $\alpha$ with the same $S$.

V-E Control strategy extraction

We conclude the section by giving a switching protocol that solves the instance (22) from a feasible solution of (24).

Theorem 3**.**

If $\{\mathbf{r}(s)\}_{s\in[T]}$, $\{\mathbf{w}(s)\}_{s\in[T+1]}$, $\{\alpha_{j}\}_{j\in[J]}$ is a feasible solution of (24), then input selection according to the switching protocol in Algorithm 1 is recursively feasible, and solves the instance (22).

Remark 3**.**

The switching protocol in Algorithm 1 depends only on the current states $\{\xi_{n}(s)\}_{n\in[N]}$ and on auxiliary information from the solution of (24). However, for implementation central coordination is required at each time step. The coordination requirement can be relaxed by simulating the system up to time $T$ and assigning an individual prefix and suffix to each subsystem. Then decentralized open-loop controllers can be constructed that realize these individual prefix-suffix paths and mimic the performance of the centralized protocol without communication requirements.

VI Analysis of the Proposed Linear Program

Theorem 3 establishes that solving (24) provides a correct solution to a deterministic discrete instance of Problem 1. We now discuss completeness of the solution approach and specify the information that can be obtained from (in)feasibility of (24) when it is solved as an integer linear program, and when integer constraints are relaxed to obtain more efficiently solvable linear programs. Table I summarizes the results of this section.

VI-A Graph Properties and Aggregate Controllability

We start by connecting properties of the induced graph $G=(\mathcal{Q},\longrightarrow)$ to a notion of reachability in the aggregate dynamics (18). These results are used in the proof of Theorem 6 below, but are also interesting in their own right since they admit a characterization of reachability in the aggregate picture.

We first define a concept of controllability on a subset of nodes $D\subset\mathcal{Q}$ for the aggregate dynamics $\Gamma$. Similarly as for controllability of linear systems on a subspace, controllability of $\Gamma$ on $D$ means that the system can be steered between any two aggregate states with support on $D$.

Definition 10**.**

A subset of nodes $D\subset\mathcal{Q}$ is completely controllable for $\Gamma$ if for any two states $\mathbf{w}$, $\mathbf{w}^{\prime}$ with support333A state $\mathbf{w}$ having support on $D$ means that $w_{q}=0$ for $q\not\in D$. on $D$ such that $\|\mathbf{w}\|_{1}=\|\mathbf{w}^{\prime}\|_{1}$, there exists a finite horizon $T$, states $\mathbf{w}(s)$, and controls $\mathbf{r}(s)$ satisfying (17), such that $\mathbf{w}(0)=\mathbf{w}$, $\mathbf{w}(T)=\mathbf{w}^{\prime}$, and $\mathbf{w}(s+1)=B\mathbf{r}(s)$ for $s\in[T]$.

Theorem 4**.**

If a strongly connected component $D$ is aperiodic, it is completely controllable for $\Gamma$ in (18).

Proof.

It is known that the incidence matrix $A_{D}$ of an aperiodic, strongly connected graph $D$ is primitive [24], i.e., there exists an integer $T$ such that all entries of $A_{D}^{T}$ are positive. This means that for each node pair $(q_{j},q_{l})$, there exists a path of length $T$ that connects them. Thus, by sending $p_{jl}$ systems along paths $q_{j}\rightarrow q_{l}$ such that $\sum_{l}p_{jl}=\mathbf{w}_{j}$ and $\sum_{j}p_{jl}=\mathbf{w}_{l}^{\prime}$, the state at time $T$ is equal to $\mathbf{w}^{\prime}$. We can define aggregate controls $\mathbf{r}(s)$ that realize these paths by switching the correct number of systems at each node over time. ∎

In the case of periodicity, it is not possible to reach every state since the parity structure of the initial state is preserved along the trajectories. However, within this restriction, the system is still controllable in a certain sense. If a strongly connected component $D$ has period $P$, its nodes can be labeled with a function $L_{P}:D\rightarrow[P]$ such that a node $q_{1}$ with $L_{P}(q_{1})=p$ only has edges to nodes $q_{2}$ with $L_{P}(q_{2})=(p+1)\mod P$. Let $D_{0},\ldots,D_{P-1}$ be the subsets of nodes induced by the equivalence relation $q_{1}\sim q_{2}$ iff $L_{P}(q_{1})=L_{P}(q_{2})$.

Corollary 1**.**

The subsets of nodes $D_{p}$ for $p\in[P]$ as constructed above are completely controllable for $\Gamma$.

Proof.

We can connect the nodes in $D_{i}$ with edges that correspond to paths of length $P$ in $D$. By construction, the resulting graphs are aperiodic, so the previous result applies. ∎

It is well known that if a discrete-time linear system is completely controllable, its reachable set for a time $s\geq n$, where $n$ is the system dimension, is the entire state-space; otherwise, it is an affine subspace that depends on $s$ and the initial state. The preceding results show a corresponding result for the aggregate dynamics $\Gamma$—a linear system with input constraints evolving on an integer lattice. The controllability in this setting is entirely characterized by the properties of the graph representing the abstraction. For the controllable case (aperiodic graph), the reachable set of $\Gamma$ at a time $s\geq T$, where $T$ is the controllability horizon from the proof of Theorem 4, is the set of all positive integer-valued vectors $\mathbf{w}$ that satisfy $\|\mathbf{w}\|_{1}=N$. In case of periodicity with a period $P$, it is the intersection of this lattice set with an affine subspace that depends on $(s\mod P)$ and the parity structure of the initial state.

Example 6**.**

Consider the graph in Figure 3, it consists of two cycles of length 2 and 4, so the period is 2. The two equivalence classes induced by the periodicity are $\{q_{0},q_{2}\}$ and $\{q_{1},q_{3},q_{4}\}$; all subsystems in the first class move to the second, and vice versa. Consider an initial state such that

[TABLE]

Due to periodicity, this parity structure is preserved in the sense that

[TABLE]

and conversely for $\mathbf{w}_{q_{1}}(s)+\mathbf{w}_{q_{3}}(s)+\mathbf{w}_{q_{4}}(s)$. However, within this restriction Corollary 1 implies that any assignment is reachable when $s$ is large enough.

VI-B Converse Results

The first result states that the restriction to prefix-suffix form is without loss of generality, provided that the prefix horizon is sufficiently large, and that the suffix cycle set is sufficiently rich.

Theorem 5**.**

Suppose that there is a solution to the instance (22). Then there is a feasible solution to (24) with a prefix length $T$ of at most $\binom{|\mathcal{Q}|+N-1}{N}$ and a suffix consisting of cycles of length at most $|\mathcal{Q}|\binom{|\mathcal{Q}|+N-1}{N}$.

The upper bounds in Theorem 5 yield large feasibility problems; next we present a result that restricts the analysis to much smaller quantities. The key observation is that an assignment can be “averaged” over its cycle without violating any counting bounds. The averaging idea is illustrated in Fig. 4 and captured in the following definition.

Definition 11**.**

For a cycle $C$, a graph period $P$ dividing $|C|$, and total weights $N_{0},\ldots,N_{P-1}$, the $P$-average assignment $\bar{\alpha}_{\{N_{p}\}_{p\in[P]}}$ is defined as

[TABLE]

In the case $P=1$, this assignment has a constant $X$-count for any cycle, more precisely;

[TABLE]

for all $s$, where $\left\langle C,\mathbf{1}\right\rangle^{X}$ simply counts the number of node-action pairs in $C$ that are in the counting set ${X}$. As a consequence, for any assignment $\alpha$, it holds that

[TABLE]

In other words, if the averaged assignment for a given total weight does not satisfy counting bounds, no assignment does.

A more general result (Lemma 1 in the appendix) shows that a cyclic integer suffix can be mapped into a (possibly non-integer) suffix defined on simple cycles via averaging. It turns out that these averaged assignments can be reached by an infinitesimal relaxation of the counting bounds since they preserve the parity structure of the initial condition. The controllability results in Section VI-A do not take counting constraints into account; when such constraints are present they may conflict with controllability. Nevertheless, by introducing an arbitrarily small relaxation of the counting constraints we can ensure that controllability is preserved. The magnitude of the relaxation does however impact the worst-case time required to control the aggregate system to a new state.

Theorem 6**.**

Suppose that there exists an integer solution to the instance $(N,\Sigma_{DFTS},\{\xi_{n}(0)\}_{n\in[N]},\{X_{l},R_{l}\}_{l\in[L]})$. Let $\operatorname{diam}(G)$ be the diameter of the induced graph. Then, if every counting constraint $(X_{l},R_{l})$ is relaxed with an absolute factor $\epsilon$ to $(X_{l},R_{l}+\epsilon)$, the non-integer version of the linear program (24) with prefix horizon $\frac{(\operatorname{diam}(G)^{2}+1)N}{\epsilon}$ and the cycle set consisting of all simple cycles, is feasible.

VI-C Rounding of Non-Integer Solution

If the linear program (24) is too large to be solvable as an integer program it may still be possible to solve it using a standard LP solver and round the result to obtain an integer solution. One option is to use a probabilistic discrepancy-minimizing rounding procedure (e.g. [25]) which allows specified relationships to be preserved after rounding; thus introducing a counting constraint violation but maintaining the validity of the solution (e.g. dynamics, prefix-suffix connection). Here we instead propose a heuristic to round only the suffix part of the solution and analyze its performance under certain assumptions on cycle structure. For a given integer suffix, the prefix part of (24) can be solved to find a matching prefix—a problem that is typically much smaller.

Counting constraints may be violated as a result of the rounding; we give bounds for the magnitude of the worst-case violation. Given an aperiodic graph and (non-integer) cycle-assignment pairs $\{C_{j},\alpha_{j}\}_{j\in[J]}$ that satisfy the counting constraints, we propose the following rounding procedure:

Step 1: Assign an integer number of subsystems to each cycle that is close to the original weight, i.e., find integers $N_{j}$ s.t. $\sum_{j\in[J]}N_{j}=\sum_{j\in[J]}\|\alpha_{j}\|_{1}$ and s.t. $N_{j}$ is close to $\|\alpha_{j}\|_{1}$. This can easily be achieved in a way s.t. $\left|N_{j}-\|\alpha_{j}\|_{1}\right|\leq 1$.

Step 2: Find individual integer assignments with total weight $N_{j}$ that are close to the average assignments, i.e., find integer assignments $\tilde{\alpha}_{j}$ s.t. $\|\tilde{\alpha}_{j}\|_{1}=N_{j}$ and s.t. $\alpha_{j}$ is close to $\bar{\alpha}_{N_{j}}$.

To this end, we let $\kappa_{1}$ and $\kappa_{2}$ be the quotient and remainder when dividing $N_{j}$ by $|C_{j}|$, i.e. $\kappa_{1}=\lfloor N_{j}/|C_{j}|\rfloor$ and $\kappa_{2}=N_{j}\mod|C_{j}|$. Then let $d=|C_{j}|/\kappa_{2}$. We consider the pseudo-periodic assignment $\tilde{\alpha}_{j}$ defined as follows:

[TABLE]

This assignment is pseudo-periodic in the sense that the 1’s are evenly distributed with distance $d$ before they are rounded to integer points, as illustrated in Fig. 5. It is easy to see that $\|\tilde{\alpha}_{j}\|_{1}=N_{j}$.

Step 3: Find a prefix using (24) that steers to $\{C_{j},\tilde{\alpha}_{j}\}_{j\in[J]}$.

If the graph is periodic the rounding can be done so as to preserve the parity structure of the original solution and guarantee reachability, but the details are omitted here. Assuming aperiodicity, we want to compare the counting bounds for the rounded solution with the counting bounds for the original solution. First we give a result that assumes a certain structure of a cycle, namely that all nodes that contribute to the counting set $X$ are placed in sequence in the cycle. Such structure is often present in practical examples, as connected counting regions tend to lead to consecutive parts in cycles.

Proposition 4**.**

Let $R_{C}^{X}=\left\langle C,\mathbf{1}\right\rangle^{X}$ be the number of nodes in a cycle $C$ that contribute to $X$-counting. If all such nodes are consecutive, i.e., $(q_{i},\mu_{i})\in{X}$ for $i\in[R_{C}^{X}]$ and $(q_{j},\mu_{j})\not\in X$ for $j\in[|C|]\setminus[R_{C}^{X}]$, then the rounding procedure (28) satisfies

[TABLE]

Proof.

When the assignment $\tilde{\alpha}_{j}$ circulates in $C_{j}$, exactly $R_{C}^{X}$ contiguous indices of $\tilde{\alpha}_{j}$ contribute to the $X$-count. We bound the number of contributing indices with value $\kappa_{1}+1$. Let $[i_{0},i_{0}+R_{C}^{X}-1]$ be a (circular) sequence representing $R_{C}^{X}$ contributing indices. Consider Fig. 5; any point that ends up in the sequence after left-rounding must satisfy $dk\in[i_{0},i_{0}+R_{C}^{X})$, where $k\in[\kappa_{2}]$. Since each left-closed, right-open interval of length $d$ captures exactly one point of the form $dk$, there are at most $\left\lceil R_{C}^{X}/d\right\rceil$ such points. Therefore by (26),

[TABLE]

∎

Corollary 2**.**

If $C_{j}$ has at most $p$ ${X}$-segments (consecutive nodes taking values in $X$), the rounding (28) satisfies $\operatorname{maxcnt}^{X}(C_{j},\tilde{\alpha}_{j})\leq\operatorname{maxcnt}^{X}(C_{j},\bar{\alpha}_{N_{j}})+p$.

Proof.

Follows from applying Proposition 4 to each of the $p$ segments. ∎

We now incorporate these results into a bound on the counting constraint violation in the overall rounding procedure. Again, this bound does not depend on the total number of subsystems $N$. The difference between the original counting bounds and their relaxed counterparts therefore becomes insignificant as $N$ grows.

Theorem 7**.**

For a given counting constraint $({X},R)$ satisfied by $\{\alpha_{j}\}_{j\in[J]}$, the following bound holds for the overall suffix rounding procedure

[TABLE]

where $p_{j}^{X}$ is the number of ${X}$-segments in the $j$’th cycle. Thus the relaxed counting constraint $({X},R+J+\sum_{j\in[J]}p_{j}^{{X}})$ is guaranteed to be satisfied by the rounded suffix $\{\tilde{\alpha}_{j}\}_{j\in[J]}$.

Proof.

By the rounding procedure, Corollary 2, and (27),

[TABLE]

Noting that $\operatorname{maxcnt}^{X}(C_{j},\bar{\alpha}_{\|\alpha_{j}\|_{1}})\leq\operatorname{maxcnt}^{X}(C_{j},\alpha_{j})$ and summing over $j\in[J]$ gives the result. ∎

If the structure required in Proposition 4 is not present, the following is a worst-case bound on the counting constraint violation due to rounding.

Theorem 8**.**

The rounding (28) satisfies

[TABLE]

Proof.

We assume that $N_{j}<|C_{j}|$, since any multiple of $|C_{j}|$ can be assigned as the average assignment. The number of nodes contributing to the $X$-count is upper bounded by $\min(R^{X}_{C},N_{j})$, hence,

[TABLE]

where the last step follows from $\max_{a,b\in[0,1]}\min(ab,(1-a)(1-b))=1/4$. ∎

These rounding bounds can be precomputed (using all cycles used in the LP instead of all cycles with non-zero assignments) and the counting constraints can be strengthened accordingly to ensure that the original constraints are satisfied after rounding.

VII Extension to Strong Heterogeneity

Up until now we considered mild heterogeneity in the continuous dynamics, which allowed us to construct a single abstraction that captures all the possible behaviors. If there is significant heterogeneity in the collection of subsystems, this may no longer be possible while maintaining a good level of approximation.

To alleviate this shortcoming, we can extend our synthesis method to a multi-class setting where each subsystem belongs to a particular class, and each class has only mild heterogeneity among its members. One abstraction per class can then be constructed, and the counting problem can be solved jointly for the different classes. In addition, counting constraints can be extended to capture class identity. For instance, we can posit that at least $R_{1}$ subsystems of class 1 be present in a given area, or guarantee that no more than $R_{2}$ subsystems of class 2 are in a particular dynamic mode.

To formalize these ideas, consider $H$ transition systems $\Sigma^{h}=(Q^{h},\mathcal{U}^{h},\longrightarrow^{h},Y^{h})$ for $h\in[H]$. Then the multi-class counting problem is as follows:

Problem 2**.**

Consider $N$ subsystems divided in $H$ classes such that class $h$ has $N_{h}$ members, and $\sum_{h\in[H]}N_{h}=N$. Subsystems in class $h$ are governed by the transition system $\Sigma^{h}$. Assume that for all $h\in[H]$, initial states $\{\mathbf{x}_{n}^{h}(0)\}_{n\in[N_{h}]}$ are given.

Given $L$ counting constraints $\{\prod_{h\in[H]}X_{l}^{h},R_{l}\}_{l\in[L]}$ with counting sets $X_{l}^{h}\subset\mathcal{Q}^{h}\times\mathcal{U}^{h}$, synthesize individual switching protocols $\{\pi_{n}^{h}\}_{n\in[N_{h}]}$ such that the generated actions $\sigma_{n}^{h}(0)\sigma_{n}^{h}(1)\sigma_{n}^{h}(2)\ldots$ and trajectories $\mathbf{x}_{n}^{h}(0)\mathbf{x}_{n}^{h}(1)\mathbf{x}_{n}^{h}(2)\ldots$ satisfy the counting constraints

[TABLE]

The linear program (24) can easily be extended to the multi-class setting, at the cost of additional variables and constraints. Assuming similar abstraction parameters and cycle selections, a problem with two classes has roughly twice as many variables as a problem with a single class. The next section includes an example showcasing how a multi-class counting problem can account for parameter heterogeneity in a family of continuous-time systems.

VIII Examples

We showcase the method on two examples; one numerical example and the TCL scheduling problem. The examples are computed with our prototype implementation available at https://github.com/pettni/mode-count, which uses Gurobi [26] as the underlying ILP solver.

VIII-A Numerical Example

Our first example is the following non-linear system:

[TABLE]

It can be shown that for a constant $\mathbf{u}$ the system is incrementally stable and that the $\mathcal{KL}$-function

[TABLE]

satisfies (7). We consider two modes $\mu_{1}$ and $\mu_{2}$ corresponding to the constant inputs $\mathbf{u}=-1$ and $\mathbf{u}=1$.

We consider the domain $\mathcal{X}=\{(x_{1},x_{2}):x_{1}\in[-2,2],x_{2}\in[-1.5,1.5]\}$. For a large $N$, we introduce mode-counting constraints $(X_{\mu_{1}},0.55N)$ and $(X_{\mu_{2}},0.55N)$ for $X_{m_{s}}=\mathcal{X}\times\{s\}$ stating that at most 55% of the subsystems can use the same dynamical mode at any given time. In addition we consider a balancing objective; namely that no more than $55\%$ of the subsystems should be in one of the sets $X_{1}=\{(x_{1},x_{2}):x_{1}\geq 0\}$ or $X_{2}=\{(x_{1},x_{2}):x_{1}\leq 0\}$; expressed by the counting constraints $(X_{1},0.55N)$ and $(X_{2},0.55N)$. Furthermore, we want all subsystems to repeatedly visit these two sets.

Using abstraction parameters $\eta=0.05$, $\tau=0.32$ we obtain an abstraction with 4,941 states that is $0.1$-approximately bisimilar to the time-discretization of (29). In order to guarantee that the counting constraints are satisfied, we therefore need to expand the counting sets as $\tilde{X}_{1}=\{(x_{1},x_{2}):x_{1}\geq-0.1\}$ and $\tilde{X}_{2}=\{(x_{1},x_{2}):x_{1}\leq 0.1\}$. We proceed by solving the discrete counting problem with randomized initial conditions and a horizon $T=10$. We sampled 200 randomized cycles that visit both $\tilde{X}_{1}^{C}$ and $\tilde{X}_{2}^{C}$, in order to achieve the second objective. We solved the problem for $N=10^{k}$ for $k=2,\ldots,9$; the solving times are shown in Table II and illustrate that the difficulty of this problem is largely independent of $N$.

Fig. 6 illustrates the number of systems that are in the counting sets over time in a trajectory, and Fig. 7 demonstrates some of the cycles that make up the suffix part of the solution.

VIII-B Application Example: TCL Scheduling

We use the following model for the dynamics of the temperature $\mathbf{x}_{n}$ of an individual TCL [15]:

[TABLE]

We assume that there are two distinct populations of TCLs, i.e., that all sets of parameters $(a_{n},b_{n},\theta_{n}^{a},P_{n}^{m})$ are $\overline{\delta}$-close to one of two nominal parameter configurations. Each nominal parameter configuration represents a mildly heterogeneous class (c.f. Section VII). The parameter values for the nominal configurations are listed in Table III along with the abstraction parameters $\eta$ and $\tau$ used for each class and the allowed deviation $\overline{\delta}$ from these nominal values. These parameters result in finite abstraction with 1,600 and 1,200 states for the two classes, respectively.

It can easily be shown that the $\mathcal{KL}$-function $\beta(r,s)=re^{-sa_{n}}$ satisfies (7) with respect to (30). In addition, it satisfies the approximate bisimulation inequality (11) for an approximation level $\epsilon=0.2$ and using the Lipschitz constant $a_{n}$ for (30); thus the results from Section IV-B apply. The constraints for the TCL problem have been introduced earlier in (4)-(5). We posit that all subsystems must remain in the temperature interval $[21.3,23.7]$; taking the approximation into account this implies that the constraint for the discrete problem becomes $\mathbf{x}_{n}(t)\in[21.5,23.5]$.

We randomly selected initial conditions for 10,000 systems of each class, sampled 50 random cycles for each class444To promote diversity in the cycle sets, the cycles were selected in order to have different fractions of time in mode on., and solved the ILP (24) for a prefix length of 20 steps (corresponding to one hour). We also introduced randomized additive model errors $\delta_{n}$ such that $|\delta_{n}|\leq\overline{\delta}$ to represent mild in-class heterogeneity. Figure 8 shows simulated trajectory densities for two different mode-on-counts, one maximal count of 6,000, i.e. $\sum_{n\in[N]}\mathds{1}_{\{\texttt{on}\}}\left(\sigma_{n}(t)\right)\in[0,6000]$; and one minimal count of 6,700, i.e. $\sum_{n\in[N]}\mathds{1}_{\{\texttt{on}\}}\left(\sigma_{n}(t)\right)\in[6700,N]$. For comparison, the fundamental minimal upper bound is 5,595 and the fundamental maximal lower bound is 7,045 as computed from formulas in [27] that apply to a centralized full-state feedback coordinator with arbitrarily fast switching. While the objective here was not to find the maximal ranges (which could be done by adding an objective function to (24)), the fundamental limits can not be attained due to approximation errors stemming from the approximate bisimulation, incomplete cycle selection, a minimal dwell time imposed by the time discretization, etc. Figure 9 shows mode-on-counts during the same simulation for the two experiments. As can be seen, the imposed counting bounds are satisfied.

IX Conclusion

This paper was concerned with control synthesis for very high-dimensional but permutation-symmetric systems subject to likewise symmetric counting constraints, and presented a scalable sound and (almost) complete solution to this problem. The main insight is to aggregate the individual subsystem dynamics as an integer linear system induced from an abstraction constructed for a single subsystem, thus avoiding abstraction of the entire state space. As we used an approximately bisimilar system as an abstraction, the same abstraction can represent not only identical subsystems (homogeneity) but also subsystems with slightly different parameters (i.e., almost symmetric, or mildly heterogeneous) by a slight change in the approximation factor. The control synthesis problem was then reduced to one of coordinating the number of subsystems that are in different parts of the discrete state-space of this abstraction. We characterized prefix-suffix solutions as the feasible set of an integer linear program, and showed how to interpret (in)feasibility of the program both in the integer and non-integer case.

The results were demonstrated on a TCL scheduling problem including tens of thousands of subsystems. With the proposed approach, it is possible to impose hard constraints on the overall power consumption of TCLs over an infinite time horizon, to the best of our knowledge this is a first in this domain. Counting constraints are also relevant in other application domains, including multi-agent planning and coordination as shown in [18].

Exploiting symmetries to achieve scalability in correct-by-construction methods is a promising direction and we will explore other types of symmetries in our future work. Another interesting direction is to consider other types of abstractions, including non-deterministic ones, since not all systems admit finite bisimulations. Although the idea of an aggregate system can still be used in this case, one should either solve a robust uncertain ILP, which could lead to conservative results, or consider reactive feedback solutions, for which different synthesis techniques should be developed.

Acknowledgment

The authors would like to thank Johanna Mathieu for insightful discussions regarding TCL coordination problem. Petter Nilsson is supported by NSF grant CNS-1239037. Necmiye Ozay is supported in part by NSF grants CNS-1446298 and ECCS-1553873, DARPA grant N66001-14-1-4045 and an Early Career Faculty grant from NASA’s Space Technology Research Grants Program.

Appendix A Proofs

Proof of Proposition 1.

For a grid point $q\in\kappa_{\eta}(\mathcal{X})$ and $x\in\mathcal{X}$, consider the relation $q\sim x$ iff $\|q-x\|_{\infty}\leq\epsilon$. Since $\epsilon>\eta/2$ this relation has the property that every $x\in\mathcal{X}$ is related to some $q\in\kappa_{\eta}(\mathcal{X})$, and conversely every $q\in\kappa_{\eta}(\mathcal{X})$ is related to at least one $x\in\mathcal{X}$. The relation satisfies 1) of Definition 2. To show that also 2) and 3) hold, we show the stronger property that if $q\sim x$, then for all $x^{\prime}=\phi_{\mu}(\tau,x,\mathbf{d})$ pertaining to an admissible $\mathbf{d}:[0,\tau]\rightarrow\mathcal{D}$ it holds that $q^{\prime}\sim x^{\prime}$ for (the unique) $q^{\prime}$ such that $q\underset{\tau,\eta}{\overset{\mu}{\longrightarrow}}q^{\prime}$.

From the construction of $\Sigma_{\tau,\eta}$, $\|q^{\prime}-\phi_{\mu}(\tau,q,\mathbf{0})\|_{\infty}\leq\eta/2$. Furthermore, under the continuity assumptions, it is known that the flow $\phi_{\mu}(t,x,\mathbf{0})$ of $f_{\mu}(x,0)$ and the flow $\phi_{\mu}(t,x,\mathbf{d})$ for any $\mathbf{d}(t)$ such that $\max_{t\in[0,\tau]}\|f_{\mu}(x,0)-f_{\mu}(x,\mathbf{d}(t))\|_{\infty}\leq\overline{\delta}_{\mu}$ satisfy $\|\phi_{\mu}(t,x,\mathbf{0})-\phi_{\mu}(t,x,\mathbf{d})\|_{\infty}\leq(\overline{\delta}_{\mu}/K_{\mu})(e^{K_{\mu}t}-1)$ [28]. Thus, it follows that for any such $\mathbf{d}:[0,\tau]\rightarrow\mathcal{D}$,

[TABLE]

Hence $q^{\prime}\sim\phi_{\mu}(\tau,x,\mathbf{d})$, which completes the proof. ∎

Proof of Theorem 1.

Let $\{\pi_{n}\}_{n\in[N]}$ be individual switching protocols that solve (12) by generating trajectories $\xi_{n}(0)\xi_{n}(1)\ldots$ and actions $\sigma_{n}(0)\sigma_{n}(1)\ldots$ for $\Sigma_{\tau,\eta}$. Due to bisimilarity and the $\eta/2$-proximity of initial conditions, the individual trajectories $\xi_{n}(0)\xi_{n}(1)\ldots$ of $\Sigma_{\tau,\eta}$ and the individual trajectories $\mathbf{x}_{n}(0)\mathbf{x}_{n}(1)\ldots$ of $\Sigma_{\tau}$ satisfy $\|\xi_{n}(s)-\mathbf{x}_{n}(s)\|_{\infty}\leq\epsilon$ for all $s\in\mathbb{N}$ when the action sequence $\sigma_{n}(0)\sigma_{n}(1)\ldots$ is implemented for both systems. By assumption,

[TABLE]

Thus for a counting set $X_{l}=X_{l}^{\mathcal{X}}\times X_{l}^{\mathcal{U}}$,

[TABLE]

where the last step follows from knowing that $\xi_{n}(s)$ only takes values $x$ such that $\kappa_{\eta}(x)=x$. Thus,

[TABLE]

which shows that the constraint in (13) is satisfied. ∎

Proof of Theorem 2.

Suppose for contradiction that there is a solution to (15) but not to (14). Let $\{\pi_{n}\}_{n\in[N]}$ be individual switching policies that solve (15) by generating trajectories $\mathbf{x}_{n}(0)\mathbf{x}_{n}(1)\ldots$ and actions $\sigma_{n}(0)\sigma_{n}(1)\ldots$ for $\Sigma_{\tau}$. Due to bisimilarity and the $\eta/2$-proximity of initial conditions, the individual trajectories $\xi_{n}(0)\xi_{n}(1)\ldots$ of $\Sigma_{\tau,\eta}$ and the individual trajectories $\mathbf{x}_{n}(0)\mathbf{x}_{n}(1)\ldots$ of $\Sigma_{\tau}$ satisfy $\|\xi_{n}(s)-\mathbf{x}_{n}(s)\|\leq\epsilon$ for all $s\in\mathbb{N}$ when the actions $\sigma_{n}(0)\sigma_{n}(1)\ldots$ are implemented for both systems. For a set $A$ we have $\kappa_{\eta}\left(A\ominus\mathcal{B}_{\infty}\left(0,\frac{\eta}{2}\right)\right)\subset\{x\in A:\kappa_{\eta}(x)=x\}\subset A$. Thus,

[TABLE]

It follows that

[TABLE]

so $\{\pi_{n}\}_{n\in[N]}$ is a solution also for (14)—a contradiction. ∎

Proof of Theorem 3.

We first consider the case $s\leq T$, and claim that the selection on Line 2 is possible if

[TABLE]

and, furthermore, that the selection guarantees that (31) holds at time $(s+1)$.

Due to (23), equation (31) holds at $s=0$. For induction, assume that (31) holds at time $s$. Then by (24e),

[TABLE]

The selection on line 2 amounts to for each $q\in\mathcal{Q}$ assigning $w_{q}(s)$ objects to $|\mathcal{U}|$ “bins” such that the $\mu$-bin has $r_{q}^{\mu}$ members; by (32) this is doable. Remark that if $\sigma_{n}(s)=\mu$, then $\xi_{n}(s+1)=q$ if and only if $\xi_{n}(s)\in\mathcal{N}_{q}^{\mu}$. Thus,

[TABLE]

where the last step follows from (16). Thus (31) holds for all $s\in[T+1]$.

Secondly, we consider the case $s\geq T$ and claim that the selection on line 4 is possible if for all $q\in\mathcal{Q}$

[TABLE]

and, furthermore, that the selection guarantees that (33) holds at time $(s+1)$. To show that (33) enables a selection $\{\sigma_{n}(s)\}_{n\in[N]}$ satisfying line 4, it suffices to remark that the selection problem is equivalent to above with $r_{q}^{\mu}(s)$ replaced by $\sum_{j\in[J]}\left\langle C_{j},\alpha_{j}^{\circlearrowleft(s-T)}\right\rangle^{\{(q,\mu)\}}$.

Due to (24c) and (31), equation (33) holds at $s=T$. Suppose for induction that (33) holds at time $s\geq T$ and that a selection $\{\sigma_{n}(s)\}_{n\in[N]}$ satisfying line 4 has been made. Then,

[TABLE]

The last step follows from the observation that a node $q^{\prime}\in\mathcal{N}_{q}^{\mu}$ must have a cycle index $(i-1)\mod|C_{j}|$ in cycle $C_{j}$, where $i$ is the cycle index of $q$. Thus the selection on line 4 is feasible for all $s\geq T$.

Finally, we show that each counting constraint $(X_{l},R_{l})$ is satisfied. For $s<T$ we have from line 2 and (24a):

[TABLE]

Thus the counting constraints are satisfied in the prefix phase. For the suffix phase, from line 4 and (24b) it follows that

[TABLE]

Thus the switching protocol in Algorithm 1 generates inputs and trajectories that satisfy the constraints of (22). ∎

Proof of Theorem 5.

Let $\{\pi^{*}_{n}\}_{n\in[N]}$ be a solution to the instance (22) and consider the generated sequence of controls $\sigma_{n}(s)$, $s\in\mathbb{N}$, and aggregate states $w_{q}(s)=\sum_{n\in[N]}\mathds{1}_{\{q\}}\left(\xi_{n}(s)\right)$. The number of possible values $\mathbf{w}(s)$ can take is finite and given by $\binom{|\mathcal{Q}|+N-1}{N}$—the number of ways in which $N$ identical objects (subsystems) can be partitioned into $|\mathcal{Q}|$ sets (nodes). We can therefore find times $T_{1},T_{2}$ with $T_{1}<T_{2}\leq\binom{|\mathcal{Q}|+N-1}{N}$ such that $\mathbf{w}(T_{1})=\mathbf{w}(T_{2})$. We show that the graph flows induced by $\{\pi^{*}_{n}\}_{n\in[N]}$ on the time interval $[T_{1},T_{2}]$ can be achieved with cycle assignments: we define a flow on a graph in a higher dimension, decompose it into flows over cycles, and project the cyclic flow onto the original graph.

Let $G=(\mathcal{Q},\longrightarrow)$ be the system graph, and define a new graph $H=(V_{H},E_{H})$. The node set $V_{H}=\underbrace{\mathcal{Q}\times\mathcal{Q}\times\ldots\times\mathcal{Q}}_{T_{2}-T_{1}\text{ times}}$ contains $T_{2}-T_{1}$ copies of each node in $\mathcal{Q}$, and copies of $q\in\mathcal{Q}$ are labeled $q_{s}$ for $s\in[T_{1},T_{2}]$. The set of edges is defined as

[TABLE]

An edge flow is induced on $H$ by $\{\pi^{*}_{n}\}_{n\in[N]}$, obtained by letting the flow along $(q_{s},\tilde{q}_{s+1})$ be the number of subsystems that traverses the edge $(q,\tilde{q})\in E$ at time $s$, and by letting the flow along $(q_{T_{2}},q_{T_{1}})$ be equal to the number of systems at $q$ at time $T_{1}$. By construction, this flow is balanced at each node (i.e. inflows equal outflows).

By the flow decomposition theorem [29, Theorem 3.5], we can then find cycles and assignments in $H$ that achieve this edge flow. As becomes evident from the proof in [29], these cycles are furthermore simple in $H$ and thus of length at most $|V_{H}|=|\mathcal{Q}|(T_{2}-T_{1})$. By projecting these cycles onto a single copy of $\mathcal{Q}$, we obtain cycles and assignments in $G$ that mimic the counting performance of $\{\pi^{*}_{n}\}_{n\in[N]}$ on the interval $[T_{1},T_{2}]$ when circulated.

We can therefore define a prefix-suffix strategy by taking as the prefix part $r_{q}^{\mu}(s)=\sum_{n\in[N]}\mathds{1}_{\{(q,\mu)\}}(\xi_{n}(s),\sigma_{n}(s)),\quad s\in[T_{1}]$, followed by a suffix part consisting of the cycles and assignments constructed as above. ∎

Lemma 1**.**

Suppose $C=C_{1}\cup C_{2}$ is a cycle that visits a node $q_{0}$ twice, so that it can be decomposed into two cycles $C_{1}=(q_{0},q_{1})\ldots(q_{{i}},q_{0})$ and $C_{2}=(q_{0},q_{{i+1}})\ldots(q_{{|C|-1}},q_{0})$.

Let $\alpha$ be an assignment to $C$ that satisfies $\operatorname{maxcnt}^{X}(C,\alpha)\leq R$, let $P$ be the graph period, and let

[TABLE]

Then the joint $X$-counts for the assignment $\bar{\alpha}_{\{N_{p}^{1}\}_{p\in[P]}}$ to $C_{1}$ and the assignment $\bar{\alpha}_{\{N_{p}^{2}\}_{p\in[P]}}$ to $C_{2}$, satisfy

[TABLE]

Proof of Lemma 1.

We have for all $i_{1}$ that

[TABLE]

and similarly for $C_{2}$. $P$ must divide both $|C_{1}|$ and $|C_{2}|$, so $((i_{1}\mod|C_{1}|)\mod P)=(i_{1}\mod P)$ for all $i_{1}$ and similarly for $|C_{2}|$. Below, for a cycle $C=(q_{0},\mu_{0},q_{1})(q_{1},\mu_{1},q_{2})\ldots$ we use the short hand notation $\mathds{1}^{C}_{X}(i)=1$ if $(q_{i},\mu_{i})\in X$ and $\mathds{1}^{C}_{X}(i)=0$ otherwise, to indicate that the $i$’th state-action pair is in the counting set $X$. For $j=1,2$ we get

[TABLE]

We have $i_{2}\mod P=(|C_{1}|+i_{2})\mod P$, and $i_{2}\mapsto|C_{1}|+i_{2}$ is a mapping from $i_{2}\in[|C_{2}|]$ to the corresponding index in $C$. We can therefore convert a sum over both $i_{1}$ and $i_{2}$ into a sum over the index $i$ of $C$.

[TABLE]

∎

Proof of Theorem 6.

By Theorem 5, we know that a correct solution must eventually lead to periodic behavior, and Lemma 1 shows that the suffix of such a solution can be mapped into a suffix on simple cycles consisting of $P$-averaged assignments, where $P$ is the period of the graph555In the case of a non-connected graph the analysis can be done separately for each strongly connected component.. What remains to show is that this latter suffix is reachable from the initial state while respecting the relaxed counting constraints.

Let $\mathbf{w}(s)$ for $s\in\mathbb{N}$ be the aggregate states for an integer solution to (22). From Lemma 1 we can obtain a set of simple cycles $I$ such that some $P$-averaged assignments to these cycles satisfy the counting bounds. In addition, these assignments have the same parity structure as $\mathbf{w}(s)$, and hence as the initial condition; therefore they are reachable from the initial condition by virtue of Corollary 1. We now propose a switching protocol to control the aggregate state to these assignments; the protocol consists in gradually steering a fraction of the systems from the original solution $\mathbf{w}(s)$ to the $P$-averaged assignments pertaining to the cycles in $I$.

We use the following notation: let $\alpha(s)$ be what remains of the correct trajectory $\mathbf{w}(s)$, let $\beta(s)$ be the fraction currently being steered towards assignments to the simple cycles, and let $\gamma(s)$ be the fraction that has already reached these assignments. We then have $\alpha(0)=\mathbf{w}(0)$, and $\beta(0)=\gamma(0)=0$. The overall system state is $\tilde{\mathbf{w}}(s)=\alpha(s)+\beta(s)+\gamma(s)$.

The protocol at time $s+1$ is as follows.

• •

If $\beta(s)=0$, set $\alpha(s+1)=\alpha(s)(1-\epsilon/\|\alpha(s)\|_{1})$, and $\beta(s+1)=\alpha(s)\epsilon/\|\alpha(s)\|_{1}$,

• •

If $\beta(s)$ has reached the average assignments to cycles in $I$, set $\beta(s+1)=0$, and $\gamma(s+1)=\beta(s)+\gamma(s)$,

• •

Otherwise, steer $\beta(s)$ toward the average assignments to cycles in $I$ .

We remark that the transitions are all properly connected and merely illustrate the transport of subsystem “weight” from $\alpha(s)$ via $\beta(s)$ to the average assignments represented by $\gamma(s)$. Transporting a mass $\epsilon$ takes at most time $(\operatorname{diam}(G)^{2}+1)$ by [30] which gives a bound $T\leq(\operatorname{diam}(G)^{2}+1)$ for the $T$ in the proof of Theorem 4. We can thus infer that the total transportation time is upper bounded by $(\operatorname{diam}(G)^{2}+1)N/\epsilon$, since an absolute weight $\epsilon$ is transported in each step.

We finally consider the counting bounds. By assumption, they are satisfied by $\mathbf{w}(s)$, and by Lemma 1 they are also satisfied once the average assignments to cycles in $I$ are reached. In the meantime, $\tilde{w}_{q}^{\mu}(s)=\alpha_{q}^{\mu}(s)+\beta_{q}^{\mu}(s)+\gamma_{q}^{\mu}(s)$. For every $s$, there is an integer $z\leq 1/\epsilon$ such that $\alpha(s)=\left(1-\frac{\epsilon z}{N}\right)\mathbf{w}(s)$, $\gamma(s)=\frac{\epsilon(z-1)}{N}\gamma_{0}(s)$, where $\gamma_{0}$ is the average assignment to the cycles. Furthermore $\|\beta_{q}^{\mu}(s)\|_{1}\leq\epsilon$ which shows that the counts are bounded as follows:

[TABLE]

We have therefore constructed a prefix-suffix solution, where the suffix part consists of simple cycles, such that the relaxed counting bounds $({X}_{l},R_{l}+\epsilon)$ are satisfied. ∎