Maximal Invariant Set Computation and Design for Markov Chains

TL;DR

This paper introduces an algorithm for computing the maximal invariant set of Markov chains with safety constraints, optimizing convergence to steady-state distributions via SDP, and demonstrates its application in decentralized swarm guidance.

Contribution

It presents a novel SDP-based method for efficient maximal invariant set computation and Markov chain synthesis with safety guarantees.

Findings

Efficient SDP formulation for invariant set computation.

Guaranteed finite determination of the invariant set.

Successful application to decentralized swarm guidance.

Abstract

We describe an algorithm for computing the maximal invariant set for a Markov chain with linear safety constraints on the distribution over states. We then propose a Markov chain synthesis method that guarantees finite determination of the maximal invariant set. Although this problem is bilinear in the general case, we are able to optimize the convergence rate to a desirable steady-state distribution over reversible Markov chains by solving a Semidefinite Program (SDP), which promotes efficient computation of the maximal invariant set. We then demonstrate this approach with a decentralized swarm guidance application subject to density upper bounds.

Equations

$$x^{+}=f(x).$$

$${\mathcal{O}}_t(f,\mathcal{X})=\{\xi \in {\mathbb{R}}^n\mid f^k(\xi )\in \mathcal{X},k=0,\dots ,t\}.$$

$${\mathcal{O}}_t(f,\mathcal{X})=\mathcal{X}\cap \mathrm{Pre}(\mathcal{X})\cap \cdots \cap {\mathrm{Pre}}^t(\mathcal{X}).$$

$$x^{+}=Mx,$$

$$\mathcal{X}=\Delta \cap \mathcal{P}(G,g),$$

$$\mathcal{O}_{t}=\mathcal{P}\left(\left[\begin{array}[]{c}G\\ GM\\ \vdots\\ GM^{t}\end{array}\right],\left[\begin{array}[]{c}g\\ g\\ \vdots\\ g\end{array}\right]\right).$$

$$\mathcal{P}\left(\left[\begin{array}[]{c}G\\ \mathbf{1}^{T}\\ -\mathbf{1}^{T}\\ -I\end{array}\right],\left[\begin{array}[]{c}g\\ 1\\ -1\\ 0\end{array}\right]\right)\subseteq\mathcal{P}(H,h).$$

$$Y_0(G-g{1}^T)M^0+\cdots +Y_t(G-g{1}^T)M^t\ge (G-g{1}^T)M^{t+1},$$

$$x^{+}=\left[\begin{array}[]{ccc}0.8&0.2&0\\ 0.2&0.2&0.9\\ 0&0.6&0.1\end{array}\right]x.$$

$${\left|\left|\left|G(M-v{1}^T{)}^{k}x\right|\right|\right|}_{\infty }\le {\left|\left|\left|G\right|\right|\right|}_{\infty }{\left|\left|\left|(M-v{1}^T{)}^k\right|\right|\right|}_{\infty }{\left|\left|\left|x\right|\right|\right|}_{\infty }$$

$$\le {\left|\left|\left|G\right|\right|\right|}_{\infty }{\left|\left|\left|(M-v{1}^T{)}^k\right|\right|\right|}_{\infty }.$$

$$M\odot \left(11^T-A_a^T\right)=0$$

$$\left[\begin{array}[]{cc}\lambda^{2}P&(M-v\mathbf{1}^{T})^{T}D^{T}\\ D(M-v\mathbf{1}^{T})&D+D^{T}-P\end{array}\right]\succeq 0.$$

$$-\lambda I⪯Q^{-1}MQ-r{r}^T⪯\lambda I,$$

$$M\ge 0,$$

$$1^{T}M=1^T,$$

$$Mv=v,$$

$$M\odot (11^T-A_a^T)=0,$$

$$-\lambda I⪯Q^{-1}MQ-r{r}^T⪯\lambda I.$$

Full text

I Introduction

I-A Background

Set theory plays an important role in robust control design [1], and the concept of invariance is critical to ensure that safety constraints are verifiably satisfied. For a given dynamical system, a set of states is positively invariant if once the system enters that set, it will never leave. Safety requirements therefore can be proven by showing that the initial condition is in a positively invariant subset of the safe region. The problem of verifying that a particular initial condition will satisfy safety constraints for all future times has been studied for linear [2][3][4] and polynomial [5][6] systems by constructing either a positively invariant set or an explicit reachable or controllable set. It is often useful to characterize the set of all such initial conditions—also called the maximal output admissible set—as the union of all positively invariant subsets of the safe region [7].

We will consider discrete-time, time-invariant Markov chains with finitely many states. The particular sequence of states typically cannot be determined in advance because of the stochastic dynamics, however the probability vector over the set of states evolves as a deterministic, linear system. There has been much research to find optimal policies for Markov decision process [8][9][10] which result in a closed-loop Markov chain. Semidefinite Programming (SDP) can aid in the design of Markov chains, specifically to constrain or optimize the mixing rate [11][12][13][14].

In this paper, we combine results from polytopic invariance analysis [15] and Markov chain synthesis with safety constraints [13]. We first specialize a set theoretic method of computing the maximal positively invariant set within a prescribed polytope to ergodic Markov chains. Under certain assumptions, the proposed algorithm requires a finite number of iterations to converge to the exact solution. We then propose an SDP-based synthesis method to compute a reversible Markov chain with sufficiently fast (or even optimal) mixing rate, subject to transition constraints, which ensures that the maximal admissible set can be exactly computed in finite time.

I-B Notation

The notation $x[k]$ is used for a time-dependent vector $x$ evaluated at time $k\in\{0,1,\dots\}$. When there is no explicit dependence on the time index, $x$ and $x^{+}$ may be used in place of $x[k]$ and $x[k+1]$, respectively. $\mathbf{1}$ is a column vector of ones, and $e_{i}$ is the $i^{\mathrm{th}}$ standard basis vector. [math] represents either the scalar [math], or a vector/matrix of all zeros. The inequality symbols $\geq$, $\leq$, $>$, and $<$ are interpreted elementwise for all entries of vectors and matrices, with $P\succ 0$ ($P\succeq 0$) indicating that matrix $P$ is positive-definite (positive-semidefinite). A probability vector $z$ satisfies $z\geq 0$ and $\mathbf{1}^{T}z=1$. The set of probability vectors of a particular dimension is called the probability simplex, and is written as $\Delta$. A Markov matrix, $M$, is a nonnegative, square matrix such that $\mathbf{1}^{T}M=\mathbf{1}^{T}$. A polyhedron is the intersection of finitely many half-spaces in $\mathbb{R}^{n}$, and is represented as $\mathcal{P}(G,g):=\{z\in\mathbb{R}^{n}\mid Gz\leq g\}$. The maximal positively invariant subset of $\mathcal{X}$ subject to $x^{+}=f(x)$ is written as $\mathcal{O}_{\infty}(f,\mathcal{X})$, or $\mathcal{O}_{\infty}(M,\mathcal{X})$ for the linear system $x^{+}=Mx$, as an abuse of notation. The binary operation $\odot$ denotes elementwise product of vectors or matrices, i.e., $(A\odot B)_{ij}=A_{ij}B_{ij}$. The spectral radius of matrix $A$, denoted $\rho(A)$, is the maximum magnitude of its eigenvalues.

II Problem Formulation

We consider a Markov chain over a finite set of states $\mathcal{S}=\{S_{1},\dots,S_{n}\}$, and transition probabilities $M_{ij}=\mathrm{Pr}(s^{+}=S_{i}\mid s=S_{j})$. This process results in the deterministic LTI system $x^{+}=Mx$ (i.e., $x[k+1]=Mx[k],\ \forall k\in\{0,1,2,\dots\}$), where $x_{i}[k]$ is the probability of being in state $S_{i}$ at time $k$, and $M$ is a Markov matrix. The transition matrix $M$ and the initial distribution $x[0]$ are treated as known, fixed quantities, so that the probability of encountering state $S_{i}$ at time $k$ is simply $e_{i}^{T}M^{k}x[0]$.

Safety constraints take the form $Gx[k]\leq g$ for all $k\in\mathbb{N}$. For example, we could specify upper and lower bounds on the probability mass of each individual state, on the sum of probabilities over a subset of states, or on the difference between probability masses of adjacent states. Having $Gx[0]\leq g$ is not sufficient to ensure that these constraints are always satisfied, since a distribution may eventually violate this constraint. Therefore it is useful to characterize the set of all initial conditions which ensure safety for all subsequent times, simplifying the analysis from checking infinitely many constraints $Gx[0]\leq g,\ Gx[1]\leq g,\dots$, to checking if $x[0]$ satisfies a finite number of inequalities.

Given a dynamical system $x^{+}=f(x)$, a set $\mathcal{Y}$ is positively invariant if $x\in\mathcal{Y}\Rightarrow f(x)\in\mathcal{Y}$. For any $\mathcal{X}\subset\mathbb{R}^{n}$, it can be shown that there exists a unique maximal positively invariant subset of $\mathcal{X}$ (sometimes simply called the maximal invariant set) $\mathcal{O}_{\infty}(f,\mathcal{X})$ containing all positively invariant subsets of $\mathcal{X}$ [16]. For a time-invariant Markov chain $x^{+}=Mx$, the problem of reducing safety constraints $Gx[k]\leq g,\ \forall k\geq 0$ to a set of conditions that depends only on $x[0]$ is equivalent to computing the maximal positively invariant subset of $\Delta\cap\mathcal{P}(G,g)$ with respect to the linear mapping $x\rightarrow Mx$. A related problem is how to design a Markov matrix on a given connected graph which ensures safety for all $x[0]\in\mathcal{X}_{0}$. This problem is more difficult to address because of the more complex relationship between the Markov matrix and the maximal invariant set, but a heuristic based on optimizing convergence rate may be obtained via semidefinite programming, as we will show in Section IV.

III Invariant Set Computation

For any asymptotically stable LTI system, there exists a positively invariant polyhedron which is bounded and nonempty [17][18]. It can also be shown that the maximal invariant subset of a polyhedron is also polyhedral, and can be computed exactly in finite time [2][19]. However, Markov chains are not asymptotically stable in the usual sense of the state approaching the origin. We first present the conceptual algorithm for computing the maximal invariant set for a general dynamical system, and then for the special case of a Markov chain with polyhedral safety constraints.

III-A Conceptual Algorithm

The problem of finding the maximal invariant set for a Markov chain is a special case of finding the maximal invariant subset of $\mathcal{X}\subset\mathbb{R}^{n}$ for the dynamical system

[TABLE]

In fact, showing that $x[0]$ is in any positively invariant subset of $\mathcal{X}$ proves safety, because all subsequent states will remain in this set. It is better still to efficiently characterize all points which lie in one of these sets, i.e., to find the maximal positively invariant subset of $\mathcal{X}$. The following lemma shows that the maximal invariant set is precisely the set of all initial conditions which remain in $\mathcal{X}$ for all time.

Lemma III.1**.**

Given the system $x^{+}=f(x)$ and set $\mathcal{X}\subset\mathbb{R}^{n}$, the safety condition $x[k]\in\mathcal{X}\ \forall k=0,1,2,\dots$ is satisfied iff $x[0]\in\mathcal{O}_{\infty}(f,\mathcal{X})$, the maximal positively invariant subset of $\mathcal{X}$.

A detailed proof of this lemma is omitted for space, but the idea is that $x[0]\in\mathcal{O}_{\infty}(f,\mathcal{X})$ implies that $\{x[0],x[1],\dots\}$ is a positively invariant subset of $\mathcal{X}$, hence $\{x[0],x[1],\dots\}$ is a subset of the maximal invariant set. Otherwise, $x[0]$ cannot be in an invariant subset of $\mathcal{X}$ because $x[k]\notin\mathcal{X}$ for some $k\geq 0$.

Since $\mathcal{O}_{\infty}(f,\mathcal{X})$ is the set of all initial conditions whose trajectories remain in $\mathcal{X}$ for all subsequent time steps, this set can be interpreted as the limit of the sequence $\mathcal{O}_{0}(f,\mathcal{X}),\ \mathcal{O}_{1}(f,\mathcal{X}),\ \dots$, where

[TABLE]

Defining the preimage of set $\mathcal{A}$ for system (1) as $\mathrm{Pre}(\mathcal{A})=\left\{\xi\in\mathbb{R}^{n}\mid f(\xi)\in\mathcal{A}\right\}$, we can express $\mathcal{O}_{t}$ as

[TABLE]

$\mathcal{O}_{\infty}(f,\mathcal{X})$ is said to be finitely determined if there exists some $t^{*}\in\mathbb{N}$ such that $\mathcal{O}_{t^{*}}(f,\mathcal{X})=\mathcal{O}_{\infty}(f,\mathcal{X})$. The following result and algorithm are slightly modified from those presented by Gilbert and Tan [2].

Lemma III.2**.**

If there is some $t^{*}$ such that $\mathcal{O}_{t^{*}}(f,\mathcal{X})\subseteq\mathcal{O}_{t^{*}+1}(f,\mathcal{X})$, then $\mathcal{O}_{\infty}(f,\mathcal{X})=\mathcal{O}_{t^{*}}(f,\mathcal{X})$.

Algorithm 1 uses the property in Lemma III.2 to compute $\mathcal{O}_{\infty}(f,\mathcal{X})$. If this algorithm terminates, then the resulting $\mathcal{O}_{t^{*}}$ is the maximal invariant set. Therefore to determine if $x$ is in $\mathcal{O}_{\infty}$, one only needs to check the finite set of conditions $x\in\mathcal{X},f(x)\in\mathcal{X},\dots,f^{t^{*}}(x)\in\mathcal{X}$. However, the condition $\mathcal{O}_{t}\subseteq\mathrm{Pre}(\mathcal{O}_{t})$ may be difficult to check for arbitrary $f$ and $\mathcal{X}$. Furthermore, this algorithm may never terminate, in which case the While loop may be stopped for some finite $t$ to give an outer approximation for $\mathcal{O}_{\infty}$ guaranteeing constraint satisfaction only for the next $t$ steps.

III-B Algorithm for Markov Chains with Polyhedral Constraints

Algorithm 1 converges for particular sets of constraints and dynamics, but there is no general guarantee that $\mathcal{O}_{\infty}$ is finitely determined, even for linear systems [2]. Fortunately, finite determination is guaranteed for ergodic Markov chains with polyhedral safety constraints. From this section onward, the set of constraints is assumed to take the form $x[k]\in\mathcal{P}(G,g)$, and $x[0]$ is required to be a probability vector. We consider the system

[TABLE]

where $M$ is a Markov matrix.

To implement Algorithm 1, we first should be able to compute $\mathrm{Pre}^{t+1}(\mathcal{X})$, to determine if $\mathcal{O}_{t}$ is a subset of $\mathrm{Pre}^{t+1}(\mathcal{X})$, and to represent $\mathcal{O}_{t+1}$ as the set intersection $\mathcal{O}_{t}\cap\mathrm{Pre}^{t+1}(\mathcal{X})$. For System (4), it is straightforward to show that $\mathrm{Pre}(\mathcal{P}(H,h))=\mathcal{P}(HA,h)$, and therefore that

[TABLE]

Since $\mathcal{O}_{t+1}=\mathcal{O}_{t}\cap\mathrm{Pre}(\mathcal{O}_{t})$, it follows by induction that if $\mathcal{O}_{t}$ is a polyhedron, then so is $\mathcal{O}_{k}$ for all $k\geq t+1$. This also demonstrates that if the constraints $\mathcal{X}$ are polyhedral, and if $\mathcal{O}_{\infty}$ is finitely determined, then $\mathcal{O}_{\infty}=\mathcal{O}_{t^{*}}$ is a polyhedron.

The following lemma demonstrates that checking the condition $\mathcal{O}_{t}\subseteq\mathrm{Pre}^{t+1}(\mathcal{X})$ in Algorithm 1 is equivalent to solving a linear feasibility problem [20][21].

Lemma III.3** (Extended Farkas’ Lemma).**

Let $\mathcal{P}_{1}=\{z\in\mathbb{R}^{n}\mid G_{1}z\leq g_{1}\}$ and $\mathcal{P}_{2}=\{z\in\mathbb{R}^{n}\mid G_{2}z\leq g_{2}\}$, and let $\mathcal{P}_{1}$ be bounded and nonempty. Then $\mathcal{P}_{1}\subseteq\mathcal{P}_{2}$ iff there exists a nonnegative matrix $Y$ such that $Yg_{1}\leq g_{2}$ and $YG_{1}=G_{2}$.

Proof.

(Sufficiency) Consider an arbitrary $x\in\mathcal{P}_{1}$, so by definition, $G_{1}x\leq g_{1}$. Both sides of this inequality can be left-multiplied by the nonnegative matrix $Y$ to obtain $YG_{1}x\leq Yg_{1}$, and thus $G_{2}x=YG_{1}x\leq Yg_{1}\leq g_{2}$. This shows that $x\in\mathcal{P}_{1}\Rightarrow x\in\mathcal{P}_{2}$, and so $\mathcal{P}_{1}\subseteq\mathcal{P}_{2}$.

(Necessity) These conditions are necessary as a direct result of strong duality. If $\mathcal{P}_{1}\subseteq\mathcal{P}_{2}$, then $p_{i}^{*}=\sup_{x\in\mathcal{P}_{1}}e_{i}^{T}G_{2}x\leq\sup_{x\in\mathcal{P}_{2}}e_{i}^{T}G_{2}x\leq e_{i}^{T}g_{2}$. Slater’s condition ensures that there is no duality gap, ensuring that the dual solutions $d_{i}^{*}$ are finite, and thus that the dual LPs are feasible. Since $-\infty<d_{i}^{*}=p_{i}^{*}\leq e_{i}^{T}g_{2}$, there must exist a minimizer $y_{i}^{*}$ in the dual domain for which $g_{1}^{T}y_{i}^{*}\leq e_{i}^{T}g_{2}$. It then follows that the matrix Y=\left[\begin{array}[]{ccc}y_{1}^{*}&\cdots&y_{m_{2}}^{*}\end{array}\right]^{T} satisfies the conditions $Yg_{1}\leq g_{2}$, $YG_{1}=G_{2}$, and $Y\geq 0$. ∎

Lemma III.3 is useful to verify that a polyhedron is positively invariant for linear dynamics, since the invariance condition $Hx\leq h\Rightarrow HAx\leq h$ is equivalent to $\mathcal{P}(H,h)\subseteq\mathcal{P}(HA,h)$. The additional structure of the Markov chain dynamics on the probability simplex can be exploited to arrive at the following results.

Theorem III.4**.**

Let $\Delta\cap\mathcal{P}(G,g)$ be nonempty. Then $\Delta\cap\mathcal{P}(G,g)\subseteq\mathcal{P}(H,h)$ iff there exists a nonnegative matrix $Y$ such that $H-h\mathbf{1}^{T}\leq Y(G-g\mathbf{1}^{T})$.

Proof.

(Sufficiency) Let $x$ be any probability vector such that $Gx\leq g$. Since $x\geq 0$, the inequality $H-h\mathbf{1}^{T}\leq Y(G-g\mathbf{1}^{T})$ can be right-multiplied by $x$ so that $(H-h\mathbf{1}^{T})x\leq Y(G-g\mathbf{1}^{T})x$, which simplifies to $Hx-h\leq Y(Gx-g)$. The right-hand side is the product of a nonnegative matrix $Y$ and a nonpositive vector $Gx-g$, so $Hx-h\leq Y(Gx-g)\leq 0$, which implies that $x\in\mathcal{P}(H,h)$.

(Necessity) The condition $\Delta\cap\mathcal{P}(G,g)\subseteq\mathcal{P}(H,h)$ is equivalent to

[TABLE]

From Lemma III.3, if $\Delta\cap\mathcal{P}(G,g)\subseteq\mathcal{P}(H,h)$, then there exists a nonnegative matrix \left[\begin{array}[]{cccc}Y&Y_{+}&Y_{-}&Y_{0}\end{array}\right] such that $YG+Y_{+}\mathbf{1}^{T}-Y_{-}\mathbf{1}^{T}-Y_{0}=H$ and $Yg+Y_{+}-Y_{-}\leq h$. The second inequality is equivalent to $(Y_{+}-Y_{-})\mathbf{1}^{T}\leq(h-Yg)\mathbf{1}^{T}$. Putting these two inequalities together and noting that $Y_{0}\geq 0$, we can verify that $YG+(h-Yg)\mathbf{1}^{T}\geq YG+Y_{+}\mathbf{1}^{T}-Y_{-}\mathbf{1}^{T}=Y_{0}+H\geq H$, which implies that $Y(G-g\mathbf{1}^{T})\geq H-h\mathbf{1}^{T}$. ∎

The next result gives a necessary and sufficient condition for invariance.

Theorem III.5**.**

Consider the system $x^{+}=Mx$ with $\mathbf{1}^{T}M=\mathbf{1}^{T}$, $M\geq 0$, and $x[0]\in\Delta\cap\mathcal{P}(G,g)$. The polyhedron $\Delta\cap\mathcal{P}(G,g)$ is positively invariant iff either $\Delta\cap\mathcal{P}(G,g)=\emptyset$, or there exists a nonnegative matrix $Y\geq 0$ such that $Y(G-g\mathbf{1}^{T})\geq(G-g\mathbf{1}^{T})M$.

Proof.

If $\Delta\cap\mathcal{P}(G,g)=\emptyset$, then no feasible trajectory exists, and the statement $x\in\emptyset\Rightarrow x^{+}\in\emptyset$ is true by the nonexistence of a counterexample.

If $\Delta\cap\mathcal{P}(G,g)$ is nonempty, then this set is positively invariant if $x\in\Delta\cap\mathcal{P}(G,g)\Rightarrow x^{+}\in\Delta\cap\mathcal{P}(G,g)$. $M$ is a Markov matrix, so $Mx$ is a probability vector for all $x\in\Delta$. Therefore all that remains to be shown is the condition $x\in\Delta\cap\mathcal{P}(G,g)\Rightarrow Mx\in\mathcal{P}(G,g)$, or equivalently $\Delta\cap\mathcal{P}(G,g)\subseteq\mathcal{P}(GM,g)$. Theorem III.4 is then applied by setting $H=GM$ and $h=g$ to obtain the necessary and sufficient condition that there exists a nonnegative matrix $Y$ such that $Y(G-g\mathbf{1}^{T})\geq GM-g\mathbf{1}^{T}$, or equivalently, $Y(G-g\mathbf{1}^{T})\geq(G-g\mathbf{1}^{T})M$.

∎

Remark 1**.**

For all $x\in\Delta$, $Gx\leq g$ is equivalent to the conical condition $(G-g\mathbf{1}^{T})x\leq 0$, therefore $g$ can be set to zero without loss of generality.

By Theorem III.4, the condition $\mathcal{O}_{t}\subseteq\mathrm{Pre}^{t+1}(\mathcal{X})$ is equivalent to the existence of nonnegative matrices $Y_{0},\dots,Y_{t}$ which satisfy the elementwise inequality

[TABLE]

which leads to Algorithm 2 for computing the maximal invariant set for a Markov chain under polyhedral safety constraints.

Figure 1 illustrates Algorithm 2 for the simple 3-dimensional example with x[k]\leq\left[\begin{array}[]{ccc}0.6&0.5&0.5\end{array}\right]^{T} and dynamics

[TABLE]

There is no guarantee that this algorithm will converge in a finite number of steps for a Markov chain with multiple eigenvalues of unit magnitude. A sufficient condition for finite determination of $\mathcal{O}_{\infty}(M,\mathcal{X})$ is that the Markov chain be ergodic (i.e., irreducible and aperiodic), which is ensured by there existing some strictly positive probability vector $v$ such that $\lim_{k\rightarrow\infty}M^{k}=v\mathbf{1}^{T}$. This property not only ensures that $v$ is an eigenvector of $M$, but also that the error $e[k]=x[k]-v$ converges to zero with convergence rate $\rho(M-v\mathbf{1}^{T})$.

Proposition III.6**.**

Let $M$ be an ergodic Markov matrix. If there exists a probability vector $v\in\Delta$ such that $Mv=v$ and $Gv<g$, then Algorithm 2 terminates after a finite number of iterations.

Proof.

It can be seen by induction that $M^{k}=(M-v\mathbf{1}^{T})^{k}+v\mathbf{1}^{T}$ for all $k\geq 1$, so the condition $GM^{k}x\leq g\ \forall k\in\mathbb{N}$ is equivalent to $Gx\leq g$ and $G(M-v\mathbf{1}^{T})^{k}x\leq g-Gv$ for all $k\geq 1$. If $Gv<g$, then there exists some scalar $\varepsilon>0$ such that $\varepsilon\mathbf{1}\leq g-Gv$. Considering the induced $\infty$-norm ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\infty}$,

[TABLE]

The last inequality follows from the fact that ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|x\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\infty}\leq 1$ for all $x\in\Delta$. $\rho(M-v\mathbf{1}^{T})<1$ since $M$ is ergodic, implying that $\lim_{k\rightarrow\infty}(M-v\mathbf{1}^{T})^{k}=0$, so there must exist a $K$ such that ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(M-v\mathbf{1}^{T})^{k}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\infty}\leq\varepsilon/{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|G\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\infty}$ for all $k\geq K$. For $k\geq K$, ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|G(M-v\mathbf{1}^{T})^{k}x\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\infty}\leq\varepsilon$, implying that $G(M-v\mathbf{1}^{T})^{k}x\leq\varepsilon\mathbf{1}\leq g-Gv$, and so $GM^{k}x\leq g$ for all $x\in\Delta$.

The stopping criterion for Algorithm 2 is equivalent to that of Algorithm 1, both terminating when $\mathcal{O}_{k}\subseteq\mathrm{Pre}^{k+1}(\mathcal{X})$. This is ensured when $k\geq K-1$ because $\mathcal{O}_{K}\subseteq\Delta\subseteq\mathrm{Pre}^{K+1}(\mathcal{X})$. ∎

This constant $K$ in the proof of Proposition III.6 is a conservative estimate for the number of times the condition $\nexists Y\geq 0\ s.t\ Y(G_{k}-g_{k}\mathbf{1}^{T})\geq GM^{k+1}-g\mathbf{1}^{T}$ must be checked in Algorithm 2. The effect of $\rho(M-v\mathbf{1}^{T})$ and $\varepsilon$ on $K$ can be found by examining the asymptotic behavior of $\|(M-v\mathbf{1}^{T})^{k}\|_{\infty}$ as $k\rightarrow\infty$. By definition, $\rho(M-v\mathbf{1}^{T})=\lim_{k\rightarrow\infty}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(M-v\mathbf{1}^{T})^{k}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\infty}^{1/k}$, and so ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(M-v\mathbf{1}^{T})^{k}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\infty}$ can be roughly approximated as $\rho(M-v\mathbf{1}^{T})^{k}$. An estimate for $K$ is obtained by solving $\rho(M-v\mathbf{1}^{T})^{K}=\varepsilon/{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|G\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\infty}$ for $K$, which results in $K\approx\frac{\log(\varepsilon/{\left|\kern-0.75346pt\left|\kern-0.75346pt\left|G\right|\kern-0.75346pt\right|\kern-0.75346pt\right|}_{\infty})}{\log(\rho)}$. The number of iterations should be lowest when $\varepsilon/{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|G\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\infty}$ is large and $\rho(M-v\mathbf{1}^{T})$ is small, indicating that all initial conditions should rapidly converge to a safe neighborhood about $v$. However, this particular approximation for $K$ tends to be very conservative in practice, often leading to much smaller $k^{*}$ than would be predicted. A more refined approach may produce a better estimate of $k^{*}$, and give greater insight into how well this algorithm can be expected to perform.

IV Constrained Markov Chain Synthesis

The previous section gave a way to compute the maximal invariant subset of a polyhedron $\Delta\cap\mathcal{P}(G,g)$ under the LTI dynamics $x^{+}=Mx$, which was finitely determined under the assumptions of Proposition III.6. We now focus on designing a Markov chain that satisfies these conditions. A Markov chain is considered to act on a directed graph with adjacency matrix $A_{a}$ defined such that $[A_{a}]_{i,j}=1$ if there is an edge from $S_{i}$ to $S_{j}$, and $[A_{a}]_{i,j}=0$ otherwise. To prevent unrealizable state transitions, the linear constraint

[TABLE]

can be enforced along with other linear relations. This equation ensures that $M_{j,i}=0$ if there is no edge from $S_{i}$ to $S_{j}$, and there is no additional constraint on $M_{j,i}$ otherwise. $M$ alternatively may be defined as a sparse matrix with a sparsity pattern that is consistent with the underlying graph. The conditions for positive invariance in Theorem III.5 are linear in $Y$ and $M$ for a given pair $(G,g)$. However, this set may be too restrictive to generate a feasible $M$ for which $\Delta\cap\mathcal{P}(G,g)$ is positively invariant. Clearly if $Mv=v$ for some $v$ in $\mathcal{P}(G,g)$, then there exists an $x[0]$ such that $GM^{k}x[0]\leq g$ for all $k\in\mathbb{N}$, namely $x[0]=v$.

As demonstrated by de Oliveira et al. [22], all eigenvalues of $M-v\mathbf{1}^{T}$ have magnitude less than or equal to $\lambda\in\mathbb{R}_{+}$ iff there exist real matrices $P\succ 0$ and $D$ satisfying the bilinear matrix inequality

[TABLE]

If $D$ is held constant, condition (11) becomes linear in $P$. Then the spectral radius $\lambda$ can then be minimized with a line search over the interval $[0,1]$ to eliminate bilinearity, allowing the spectral radius to be minimized over Markov chains by solving a sequence of linear matrix inequalities (LMIs). However in the special case that $M$ is reversible, i.e., $M\mathrm{diag}(v)=\mathrm{diag}(v)M^{T}$, this condition may be replaced by the LMI

[TABLE]

with $r=v^{1/2}$ elementwise and $Q=\mathrm{diag}(r)$ [13]. Although this condition is necessary and sufficient for reversible Markov matrices, it is not necessary in the general case, and may fail to find a feasible solution if one exists.

IV-A Synthesis Procedure

Let the following information be given: a polytope of constraints $\mathcal{P}(G,g)$, and a graph with adjacency matrix $A_{a}$ which satisfies the following assumption.

Assumption 1**.**

There is some positive integer $l$ for which there exists a path of length $l$ between any two nodes.

Assumption 1 is necessary for a primitive $M$ to exist, ensuring ergodicity. If the steady-state distribution is specified as some particular $v>0$, then the constraint $Mv=v$ is linear. If $v$ is not specified, then a second assumption is made.

Assumption 2**.**

There exists some $v\in\Delta$ such that $v>0$ and $Gv<g$.

If such a $v$ exists, then $D$ can be set to $\mathrm{diag}(v)^{-1}$ as before. These two assumptions ensure that there exists some Markov matrix $M$ for which $\lim_{k\rightarrow\infty}M^{k}=v\mathbf{1}^{T}$, and that Algorithm 2 terminates as a result of Proposition III.6 from the existence of some $\varepsilon>0$ such that $g-Gv\geq\varepsilon\mathbf{1}$. Furthermore if $A_{a}$ is symmetric, then an ergodic, reversible Markov chain can be explicitly constructed, e.g., with the Metropolis-Hastings algorithm, which would ensure that a feasible solution exists [12].

If the steady-state distribution is not given, then the first step is to choose a $v\in\Delta$ that satisfies the strict inequalities $v>0$ and $Gv<g$. Then find the smallest $\lambda\geq 0$ for which there exists $M\in\mathbb{R}^{n\times n}$ subject to:

[TABLE]

If no reversible ergodic Markov matrix is found, then this problem can be modified by replacing (13e) with LMI (11) and $P\geq 0$. Bilinearity can be avoided by fixing $D$, e.g., $D=\mathrm{diag}(v)^{-1}$ as suggested in [12], and testing various values of $\lambda$ in the interval $[0,1]$ to find the minimum.

This procedure requires solving a sequence of linear feasibility problems in the matrix variable $M$. Convex programming solvers can minimize a convex function of $M$ subject to constraints (13), and additional convex constraints can be imposed such as upper or lower bounds on each element of $M$ to directly specify a range of desirable transition probabilities. One possibility is to minimize the frequency of state transitions, $\sum_{i=1}^{n}(1-M_{i,i})v_{i}$, to discourage the system from changing states too often. This quantity is the probability of a transition occurring while the system has converged to $v$.

Directly solving for an optimal $v$ may not scale well with dimension, as it appears bilinearly in the constraints. There may be a tradeoff in performance between not letting any element of $v$ from being too close to [math], and maximizing the smallest entry of $g-Gv$. The maximal invariant set computation may not terminate if $M$ has multiple eigenvalues of magnitude $1$. We minimize $\rho(M-v\mathbf{1}^{T})$ as a tractable heuristic to compute $\mathcal{O}_{\infty}$ in as few steps as possible, however further research may better estimate the number of steps the algorithm will terminate given particular $M$ and $\mathcal{P}(G,g)$.

V Example

This example is a decentralized swarm guidance problem, where there are many agents which each have the same stochastic control scheme based on a state-dependent probability distribution. Each agent must choose which control action to perform based only on knowledge of its current location, and limited collision avoidance capabilities. The sequence of bins each agent visits follows a Markov chain, and if all agents follow the same randomized policy, $x[k]$ can be interpreted as the expected distribution of agents at step $k$. The Markov chain must be designed to converge safely to a desired steady-state distribution.

The bins are arranged in a grid, with some bins containing obstacles which must be avoided. No cell may have more than 30% probability mass at any time. The steady-state distribution is such that each agent is expected to spent 90% of its time in one of the four terminal states, regardless of their initial location. The steady-state distribution $v$ is $22.5\%$ at each of the four terminal bins and $0.25\%$ elsewhere. The convergence rate was optimized by checking the feasibility of the set defined in (13) using SeDuMi, with $\lambda^{*}=0.9950$ being optimal. With this value of $\lambda$, the particular $M$ was found by maximizing the linear objective function $\sum_{i}M_{i,i}v_{i}$, the rate at which agents remain in their same bin. Algorithm 2 was run to determine set of safe initial conditions, terminating with $t^{*}=3$. We verify that the distribution in Figure 2 is in $\mathcal{O}_{\infty}$ by checking that $x[k]\leq(0.3)\mathbf{1},\ \forall k\in\{0,1,2,3\}$. This distribution is propagated forward to $k=100$ and $k=1000$, as shown in Figure 3, which confirm that the safety constraints are always satisfied.

VI Conclusion

We have specialized a general maximal invariant set computation algorithm for Markov chains subject to polyhedral safety constraints on the probability distribution over states, and we have given conditions which ensure that the maximal invariant set is finitely determined and polyhedral. We then gave an SDP-based procedure for synthesizing a Markov chain which guarantees finite determination by proper choice of the steady-state distribution and by minimizing the second largest eigenvalue of the Markov matrix $M$, by promoting rapid convergence to a safe steady-state distribution. This method was illustrated with a swarm exploration example.