{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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.\n",
"\n",
"We consider the example given in equation (15) of [RKKM05]:\n",
"$$x^+ =\n",
"\\begin{bmatrix}\n",
" 1 & 1\\\\\n",
" 0 & 1\n",
"\\end{bmatrix}x +\n",
"\\begin{bmatrix}\n",
" 1\\\\\n",
" 1\n",
"\\end{bmatrix} u\n",
"+ w$$\n",
"with the state feedback control $u(x) = -\\begin{bmatrix}1.17 & 1.03\\end{bmatrix} x$.\n",
"The controlled system is therefore\n",
"$$x^+ =\n",
"\\left(\\begin{bmatrix}\n",
" 1 & 1\\\\\n",
" 0 & 1\n",
"\\end{bmatrix} -\n",
"\\begin{bmatrix}\n",
" 1\\\\\n",
" 1\n",
"\\end{bmatrix}\n",
"\\begin{bmatrix}1.17 & 1.03\\end{bmatrix}\\right)x\n",
"+ w =\n",
"\\begin{bmatrix}\n",
" -0.17 & -0.03\\\\\n",
" -1.17 & -0.03\n",
"\\end{bmatrix}x\n",
"+ w.$$\n",
"\n",
"[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."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2×2 Array{Float64,2}:\n",
" -0.17 -0.03\n",
" -1.17 -0.03"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A = [1 1; 0 1] - [1; 1] * [1.17 1.03]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The set of disturbance is the unit ball of the infinity norm."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"V-representation Polyhedra.PointsHull{Float64,Array{Float64,1},Int64}:\n",
"4-element iterator of Array{Float64,1}:\n",
" [-1.0, -1.0]\n",
" [-1.0, 1.0]\n",
" [1.0, -1.0]\n",
" [1.0, 1.0]"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using Polyhedra\n",
"Wv = vrep([[x, y] for x in [-1.0, 1.0] for y in [-1.0, 1.0]])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We will use the default library for this example but feel free to pick any other library from [this list of available libraries](https://juliapolyhedra.github.io/) such as [CDDLib](https://github.com/JuliaPolyhedra/CDDLib.jl).\n",
"The LP solver used to detect redundant points in the V-representation is [GLPK](https://github.com/JuliaOpt/GLPK.jl). Again, you can replace it with any other solver listed [here](http://www.juliaopt.org/JuMP.jl/dev/installation/#Getting-Solvers-1) that supports LP."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"DefaultLibrary{Float64}(OptimizerFactory(GLPK.Optimizer, (), Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}()))"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using GLPK\n",
"using JuMP\n",
"lib = DefaultLibrary{Float64}(with_optimizer(GLPK.Optimizer))"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Polyhedron DefaultPolyhedron{Float64,Polyhedra.Intersection{Float64,Array{Float64,1},Int64},Polyhedra.Hull{Float64,Array{Float64,1},Int64}}:\n",
"4-element iterator of Array{Float64,1}:\n",
" [-1.0, -1.0]\n",
" [-1.0, 1.0]\n",
" [1.0, -1.0]\n",
" [1.0, 1.0]"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"W = polyhedron(Wv, lib)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The $F_s$ function of equation (2) of [RKKM05] can be implemented as follows."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Fs (generic function with 2 methods)"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function Fs(s::Integer, verbose=1)\n",
" @assert s ≥ 1\n",
" F = W\n",
" A_W = W\n",
" for i in 1:(s-1)\n",
" A_W = A * A_W\n",
" F += A_W\n",
" if verbose ≥ 1\n",
" println(\"Number of points after adding A^$i * W: \", npoints(F))\n",
" end\n",
" removevredundancy!(F)\n",
" if verbose ≥ 1\n",
" println(\"Number of points after removing redundant ones: \", npoints(F))\n",
" end\n",
" end\n",
" return F\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Number of points after adding A^1 * W: 16\n",
"Number of points after removing redundant ones: 8\n",
"Number of points after adding A^2 * W: 32\n",
"Number of points after removing redundant ones: 12\n",
"Number of points after adding A^3 * W: 48\n",
"Number of points after removing redundant ones: 16\n",
" 21.374037 seconds (60.75 M allocations: 3.034 GiB, 7.97% gc time)\n"
]
},
{
"data": {
"text/plain": [
"Polyhedron DefaultPolyhedron{Float64,Polyhedra.Intersection{Float64,Array{Float64,1},Int64},Polyhedra.Hull{Float64,Array{Float64,1},Int64}}:\n",
"16-element iterator of Array{Float64,1}:\n",
" [-1.29, -2.56]\n",
" [-0.71, 2.56]\n",
" [-0.95, 1.78]\n",
" [-0.898, 2.032]\n",
" [-0.9022, 2.0158]\n",
" [-0.938, 1.852]\n",
" [-0.77, 2.5]\n",
" [-1.29, -0.56]\n",
" [1.29, 0.56]\n",
" [0.77, -2.5]\n",
" [0.938, -1.852]\n",
" [0.9022, -2.0158]\n",
" [0.898, -2.032]\n",
" [0.95, -1.78]\n",
" [0.71, -2.56]\n",
" [1.29, 2.56]"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"@time Fs(4)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The Figure 1 of [RKKM05] can be reproduced as follows:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using Plots\n",
"plot()\n",
"for i in 10:-1:1\n",
" plot!(Fs(i, 0))\n",
"end\n",
"# The cell needs to return the plot for it to be displayed\n",
"# but the `for` loop returns `nothing` so we add this dummy `plot!` that returns the plot\n",
"plot!()"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using Plots\n",
"plot()\n",
"for i in 10:-1:10\n",
" plot!(Fs(i, 0))\n",
"end\n",
"# The cell needs to return the plot for it to be displayed\n",
"# but the `for` loop returns `nothing` so we add this dummy `plot!` that returns the plot\n",
"plot!()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now, suppose we want to compute an invariant set by scaling $F_s$ by the appropriate $\\alpha$.\n",
"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$."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"αo (generic function with 1 method)"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function αo(s)\n",
" A_W = A^s \\ W # This triggers the computation of the H-representation of W\n",
" # A_W is H-represented\n",
" hashyperplanes(A_W) && error(\"HyperPlanes not supported\")\n",
" return maximum([Polyhedra.support_function(h.a, W) / h.β for h in halfspaces(A_W)])\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We obtain $\\alpha \\approx 1.9 \\cdot 10^{-5}$ like in [RKKM05]."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.919070000000001e-5"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"α = αo(10)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The scaled set is is the following:"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using Plots\n",
"plot((1 - α)^(-1) * Fs(10, 0))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.1.0",
"language": "julia",
"name": "julia-1.1"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.1.1"
}
},
"nbformat": 4,
"nbformat_minor": 2
}