In this example, we implement the method presented in [RKKM05] to compute a robust positively invariant polytope of a linear system under a disturbance bounded by a polytopic set.

We consider the example given in equation (15) of [RKKM05]: $$ x^+ = \begin{bmatrix} 1 & 1\\ 0 & 1 \end{bmatrix}x + \begin{bmatrix} 1\\ 1 \end{bmatrix} u

- w $$ with the state feedback control $u(x) = -\begin{bmatrix}1.17 & 1.03\end{bmatrix} x$. The controlled system is therefore $$ x^+ = \left(\begin{bmatrix} 1 & 1\\ 0 & 1 \end{bmatrix} - \begin{bmatrix} 1\\ 1 \end{bmatrix} \begin{bmatrix}1.17 & 1.03\end{bmatrix}\right)x
- w = \begin{bmatrix} -0.17 & -0.03\\ -1.17 & -0.03 \end{bmatrix}x
- w. $$

[RKKM05] Sasa V. Rakovic, Eric C. Kerrigan, Konstantinos I. Kouramas, David Q. Mayne *Invariant approximations of the minimal robust positively Invariant set*. IEEE Transactions on Automatic Control 50 (**2005**): 406-410.

In [1]:

```
A = [1 1; 0 1] - [1; 1] * [1.17 1.03]
```

Out[1]:

The set of disturbance is the unit ball of the infinity norm.

In [2]:

```
using Polyhedra
Wv = vrep([[x, y] for x in [-1.0, 1.0] for y in [-1.0, 1.0]])
```

Out[2]:

We will use the default library for this example but feel free to pick any other library from this list of available libraries such as CDDLib. The LP solver used to detect redundant points in the V-representation is GLPK. Again, you can replace it with any other solver listed here that supports LP.

In [3]:

```
using GLPK
using JuMP
lib = DefaultLibrary{Float64}(GLPK.Optimizer)
W = polyhedron(Wv, lib)
```

Out[3]:

The $F_s$ function of equation (2) of [RKKM05] can be implemented as follows.

In [4]:

```
function Fs(s::Integer, verbose=1)
@assert s ≥ 1
F = W
A_W = W
for i in 1:(s-1)
A_W = A * A_W
F += A_W
if verbose ≥ 1
println("Number of points after adding A^$i * W: ", npoints(F))
end
removevredundancy!(F)
if verbose ≥ 1
println("Number of points after removing redundant ones: ", npoints(F))
end
end
return F
end
```

Out[4]:

We can see below that only the V-representation is computed. In fact, no H-representation was ever computed during `Fs`

. Computing $AW$ is done by multiplying all the points by $A$ and doing the Minkowski sum is done by summing each pair of points. The redundancy removal is carried out by CDD's internal LP solver.

In [5]:

```
@time Fs(4)
```

Out[5]:

The Figure 1 of [RKKM05] can be reproduced as follows:

In [6]:

```
using Plots
plot()
for i in 10:-1:1
plot!(Fs(i, 0))
end
```

The cell needs to return the plot for it to be displayed
but the `for`

loop returns `nothing`

so we add this dummy `plot!`

that returns the plot

In [7]:

```
plot!()
```

Out[7]:

Now, suppose we want to compute an invariant set by scaling $F_s$ by the appropriate $\alpha$.
In equation (11) of [RKKM05], we want to check whether $A^s W \subseteq \alpha W$ which is equivalent to $W \subseteq \alpha A^{-s} W$.
Note that `A^s \ W`

triggers the computation of the H-representation of `W`

and `A_W`

is H-represented.

In [8]:

```
function αo(s)
A_W = A^s \ W
hashyperplanes(A_W) && error("HyperPlanes not supported")
return maximum([Polyhedra.support_function(h.a, W) / h.β for h in halfspaces(A_W)])
end
```

Out[8]:

We obtain $\alpha \approx 1.9 \cdot 10^{-5}$ like in [RKKM05].

In [9]:

```
α = αo(10)
```

Out[9]:

The scaled set is is the following:

In [10]:

```
using Plots
plot((1 - α)^(-1) * Fs(10, 0))
```

Out[10]:

*This notebook was generated using Literate.jl.*