{
"cells": [
{
"cell_type": "markdown",
"source": [
"# Minimal Robust Positively Invariant Set"
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"### Introduction"
],
"metadata": {}
},
{
"cell_type": "markdown",
"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",
"$$\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",
"$$\n",
"with the state feedback control $u(x) = -\\begin{bmatrix}1.17 & 1.03\\end{bmatrix} x$.\n",
"The controlled system is therefore\n",
"$$\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",
"\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."
],
"metadata": {}
},
{
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": "2×2 Matrix{Float64}:\n -0.17 -0.03\n -1.17 -0.03"
},
"metadata": {},
"execution_count": 1
}
],
"cell_type": "code",
"source": [
"A = [1 1; 0 1] - [1; 1] * [1.17 1.03]"
],
"metadata": {},
"execution_count": 1
},
{
"cell_type": "markdown",
"source": [
"The set of disturbance is the unit ball of the infinity norm."
],
"metadata": {}
},
{
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": "V-representation Polyhedra.PointsHull{Float64, Vector{Float64}, Int64}:\n4-element iterator of Vector{Float64}:\n [-1.0, -1.0]\n [-1.0, 1.0]\n [1.0, -1.0]\n [1.0, 1.0]"
},
"metadata": {},
"execution_count": 2
}
],
"cell_type": "code",
"source": [
"using Polyhedra\n",
"Wv = vrep([[x, y] for x in [-1.0, 1.0] for y in [-1.0, 1.0]])"
],
"metadata": {},
"execution_count": 2
},
{
"cell_type": "markdown",
"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://jump.dev/JuMP.jl/dev/installation/#Getting-Solvers-1) that supports LP."
],
"metadata": {}
},
{
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": "Polyhedron DefaultPolyhedron{Float64, Polyhedra.Intersection{Float64, Vector{Float64}, Int64}, Polyhedra.Hull{Float64, Vector{Float64}, Int64}}:\n4-element iterator of Vector{Float64}:\n [-1.0, -1.0]\n [-1.0, 1.0]\n [1.0, -1.0]\n [1.0, 1.0]"
},
"metadata": {},
"execution_count": 3
}
],
"cell_type": "code",
"source": [
"using GLPK\n",
"using JuMP\n",
"lib = DefaultLibrary{Float64}(GLPK.Optimizer)\n",
"W = polyhedron(Wv, lib)"
],
"metadata": {},
"execution_count": 3
},
{
"cell_type": "markdown",
"source": [
"The $F_s$ function of equation (2) of [RKKM05] can be implemented as follows."
],
"metadata": {}
},
{
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": "Fs (generic function with 2 methods)"
},
"metadata": {},
"execution_count": 4
}
],
"cell_type": "code",
"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"
],
"metadata": {},
"execution_count": 4
},
{
"cell_type": "markdown",
"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."
],
"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",
" 1.575101 seconds (3.98 M allocations: 240.741 MiB, 4.29% gc time, 99.56% compilation time)\n"
]
},
{
"output_type": "execute_result",
"data": {
"text/plain": "Polyhedron DefaultPolyhedron{Float64, Polyhedra.Intersection{Float64, Vector{Float64}, Int64}, Polyhedra.Hull{Float64, Vector{Float64}, Int64}}:\n16-element iterator of Vector{Float64}:\n [-1.29, -2.56]\n [-1.29, -0.5599999999999999]\n [-0.9500000000000002, 1.7799999999999996]\n [-0.9380000000000002, 1.8519999999999996]\n [-0.9022000000000001, 2.0157999999999996]\n [-0.8980000000000001, 2.0319999999999996]\n [-0.77, 2.4999999999999996]\n [-0.71, 2.56]\n [1.29, 2.56]\n [1.29, 0.5599999999999999]\n [0.9500000000000002, -1.7799999999999996]\n [0.9380000000000002, -1.8519999999999996]\n [0.9022000000000001, -2.0157999999999996]\n [0.8980000000000001, -2.0319999999999996]\n [0.77, -2.4999999999999996]\n [0.71, -2.56]"
},
"metadata": {},
"execution_count": 5
}
],
"cell_type": "code",
"source": [
"@time Fs(4)"
],
"metadata": {},
"execution_count": 5
},
{
"cell_type": "markdown",
"source": [
"The Figure 1 of [RKKM05] can be reproduced as follows:"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"using Plots\n",
"plot()\n",
"for i in 10:-1:1\n",
" plot!(Fs(i, 0))\n",
"end"
],
"metadata": {},
"execution_count": 6
},
{
"cell_type": "markdown",
"source": [
"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"
],
"metadata": {}
},
{
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": "Plot{Plots.GRBackend() n=10}",
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO3dd2BUZb4+8Pec6ZmZdJKQQkKVQBqEkNAD0iE0payuIrrKuqvuz3X3Wnavd/eubrG7lrt6VcSKgqD0FkJPIYVQAgkJJb2RhEwv55zfH9mbmQktkMm8U57PX5nh5MzDEeeZd+Z85zCCIBAAAABfxdIOAAAAQBOKEAAAfBqKEAAAfBqKEAAAfBqKEAAAfBqKEAAAfBqKEAAAfBqKEAAAfBqKEAAAfBqKEAAAfJpLi/DLL78sKSlx5SO6nsVioR3BY+BY9R6OVe/hWPWe1WqlHcEtuLQIs7OzT5065cpHdD2j0Ug7gsfAseo9HKvew7HqPaPRiK+bJnhrFAAAfByKEAAAfBqKEAAAfBqKEAAAfBqKEAAAfBqKEAAAfBqKEAAAfBqKEAAAfJqYdgCgqbKyctmyZZcvX6by6IIgMAxD5aE9Do5V7+FY9Z4nHiupVFpbWyuXy524TxShLyotLV25cuWVSxc5geNlChGt9wUYQgi+1aJ3cKx6D8eq9zzwWOmvdeTm5k6fPt2J+0QR+pCCgoKHH374yqVLVsHKS2RELBKNHc+Ny+ACAmhHAwDoheIidu8up+8VRej9Dh8+/Nhjj9XWXOF4jpfJBTHLjExmlGq/S1eMBYXyopO0AwIA3J7FbOIyM/tjzyhCb7Zs2bKdO7ZbBZ7I5YxMRuJHEbFMVV3LVl6eEzl68Yi5M6ePEjM4YQoAPMDK3PX7+udyGShCr/W73/3ux127REo/ITGZFUllZWWyc5WLo5OWj148fsAQ1tM+IQcA6CcoQu+Uk5Pz1vvviQfFia+2q8+UL49OuW/cz1NCYmjnAgBwOyhCL9TZ2Tln7txgsXy139ClI5NHB0XRTgQA4L5QhF4oJDTEj5Hsm/NsnDq0xx9xAn+uo8HK81SCAQDcEYYho4Oi+vtUBhSht/EP8JcS0cZ7n7RvQQvPHWms+LG6cGfNqUA/sVwiopgQAKCXajsM76U/lBWb0q+PgiL0KikpKZzB/EbGyrQBgwkhJs6a11y1s+7k5ktFgUpRZoLsf+ZFxIRIaMcEAOiVP21oswr9/g4WitB7PPnkkxfKyh+Lz1w6KGV37Zkt1QV7a8oGD5BlJso+mx8xwB//rQEAbgBPjl5iy5Ytn3/2+bCAsEuddcM3Pj86SjV1tOTrZVHBKrwLCgBwKyhCb9DY2Lhy1c8CJRKLqC0xxfTrVTEBfug/AIBeQRF6g5hBg+SMSCrm3n5kYKASFQgAcAfw9Voez0/pJxaIWMS/9lAYWhAA4E5hRejZhg8fLpisCin7x/tCBodJaccBAPA8KEIPtnr16rortYEK6c+nqSaMUNKOAwDgkVCEnurbb7/duOH7ILl0crxkWbo/7TgAAJ4KnxF6pLq6utWPrFGJxYMGME/PD6IdBwDAg6EIPVLs4MEywijl/H+vChWxuKASAMDdQxF6HrlcLhEYmVj4x0NhKjn+CwIA9AmeRj1MbGwsY+UVEubPKwfgW0MBAPoOJ8t4kqVLl7Y0NKnlkkdnqMcOVtCOAwDgDVCEHuPDDz/cs3N3kFw6PVG2OA2niQIAOAfeGvUMFRUVv332twFSSVw4++RsnCYKAOA0WBF6AKvVOjohwY8VK+Tcn1aEsY6vXixW4cRFfXGN2Srg9FEA8AAMEValqcMD3KWA3CUH3IJSpZQIjETEv/7QwO7TRE0WIa9St+ucMb9CK40O0yWn8lL81wQADyA9VJAcaQwPUNEO8m946nR3AwdGsFZBJmX/9mBYZJDEZBFOXNTvKDPmndeIo8P0k6YKT6aYgvCRIQB4DElFFSHXaKewQRG6tblz51672q6WSx6f6d+h5178sS3vXKc4Okw/aYrw5BhToJp2QAAAj4cidF9vvvnmkZyDATIpUZHX91yVjBqsnzCNPD3apMTgBACA06AI3dSZM2deeukPITKp1I+rHzrc8sRKi0JGOxQAgBfC+IQ7slgsY1JTFSyrVvCtes5y/1yCFgQA6B8oQnekVCqlApGJhXtTFMLwWBIVRjsRAIDXQhG6ndDQUJYnMjHz6gMDNpca9FkzaScCAPBmTviMUBCEgoKC0tJShUIxc+bMgQMH9n2fPmvKlCn6To1KJv7NgsBOPdcp9SOjh9IOBQDgzZywInzxxRfXrFlTVFS0Y8eOkSNHHjx4sO/79E3r168vKSgMksvuy1DNTFR9lq8zLsFyEACgfzlhRfjMM8/87W9/YxiGEPLSSy/9/e9/z8zM7PtufdCTT6wNlMniY9jV0wIvNpurWqzChDG0QwEAeDknFGFkZGT3z0FBQYIg9H2fPsjPTyFj2UAV/4f7whiGrDuuscybRsQi2rkAALycM+cIW1pa3nnnnffff/9mGzQ1NW3fvr2urq7rpkwm+/Wvfy0SedVzvcVisVgsd/pby5Yt481WqUL09wcj5BKmXccdLddyv87oj4QAAB7tjp5mRSIRy97mQ0CnFaFWq12yZMny5cuXLl16s204jtPpdO3t7V035XI5x3Fd76l6DZ7neZ6/o1+5fPly9p59cgnzxsPhA/zFhJBv8zuFyWOJWtk/GQEAPJggCL1/mr1tCxJnFaFer8/Kyho9evTbb799i80iIyOnT5++evVqpzyoezKbzTLZnQ2/JyYmSlnhpaVhwyJkhBCTRdhcpDH/ZVr/BAQA8GxSqfROn2ZvzQlnjZpMpuXLl8fGxv7rX//ysuWdC0ydOlXC8GumBU+J//f6b+fJTjIijkRiiB4AwBWcUIQvvfTS3r17NRrNqlWrVqxY8dRTT/V9nz7iwoULJQW50+LVD0wO6LpHEMj6Ap1+EaYmAABcxAlvjf7sZz/LyLCd1qFU4pOt3ho1apRKSp5fHNp9z/EKnVahIvFDKKYCAPApTijCcePGjRs3ru/78TVTp06Vi4Qn7g2Vim3vJ6/L1+oXZVFMBQDga/Bdo3RUVFScyMtXSEULxtourlvVZLp4lSMTkikGAwDwNShCOhISEvwk7BMzgyUiu+XgcY1lfibxrsFKAAA3hyKkICMjQ8QTPykzN1nVfWerxnq0XMfdiyF6AACXwhXqXa2iouJMyanhQSFzMngRa1sOflegIVNTicqPYjYAAB+EFaGrjU5MUIpldbr2eSm2TweNFmFLsca8YDrFYAAAvgkrQpfKzMyUEVFSSET04FaF1PYqZEdJJxk5mISHUMwGAOCbsCJ0qSO5uUEK9Ynmi0vG25aDgkC+KNDhSvQAAFRgReg68fHxfkQ0KTROq6oK87cd+aPlOp1STUYOppgNAMBnYUXoIjzPV1y8OEAVeKT53IpJKvs/WpeP71QDAKAGRegiw4YN82PEMwYMiw4RDY+wfW96ZZPpcpuVZCRRzAYA4MtQhK7AcdyVhnqxSFzYXrlissN3sX52TGNeOB1D9AAAtKAIXSE1NVXGM4tikjqs1zKG24qwVWM9fkHPz8AQPQAANShCVzh7vkwgpF7fvHKS2m6Gnnyb30mmphKlgl40AABfhyLsd7/85S8Jz8wblFTUemWO3XeqGS3Cj8Va8/xMetEAAABF2P/WrftMwrASgVs6Xm0/RL+tuJOMGoIhegAAulCE/WvDhg1WgYyLGLan7sxSuyF6XiBfntAZsu6lmA0AAAgG6vvbI4884ieSDlX4B45UhajthujPa/VKNbkHQ/QAAJRhRdiPqqqqzIIQoQrcU3dq+UTHqYl8nX7xLFrBAACgG4qwH6WlpalY8bSQIbFh4mF2Q/Tn603V1wSSjiF6AAD6UIT9heO4awadmBXnXS1fPsnhKoOf52nM86cREQ4+AAB9eC7uLykpKTKeXRSTpON16cNs74u2dFrzLuj5eydQzAYAAN1QhP3l/IVyhpAaXdOqKSr7Ifpv8juFaeOIn5xeNAAAsEER9otf/OIXDM/Mi00pbauek2SbmjCY+a0lWsu8aRSzAQCAPRRhv/jiy/VihmV485LxapnEth7cWqwho4diiB4AwH2gCJ3viy++4ARm/MDh++rOLhnv330/L5AvC7QYogcAcCsYqHe+J5543E8kjZOrA+KVISrb9ZUOn9Ma1P5kRBy9aAAA0BNWhE5WWVlpJmSgOmhP3an7JzpeiT5PZ1gym1YwAAC4IRShk02fPl3JiKcGDx4cJh4WbjdEX2es0QpkfCLFbAAAcD0UoTNZrVaN0SATS3Nby3teiT5Pa54/HUP0AADuBs/LzpSUlCQX2IVRCSZGN36o7dtkmjutJyr1/Ix0itkAAOCGUITOdOFiJSFMta5xxWQVYz9En9fJTx+PIXoAADeEInSaNWvWsDyzIDb5VFuN/RC93sRvLdFY5k6lmA0AAG4GReg0X3/zlYhhec64dLy/VOwwRM8kDidhwRSzAQDAzaAInePTTz/lBDZ94Ij99WVL7K9Ez5OvTmCIHgDAfWGg3jl+/etfKVhxrEwVOEodbDdEf+ic1hgURIbHUswGAAC3gBWhE5w/f95MyKCAAbvrSpdPchii/yxfp180k1YwAAC4LRShE2RkZCgZ8cSgQcMiJIMHSLvvP1dnrNMKJC2BYjYAALg1FGFfWSwWjckoF0uPt55f7jhE/2mu1rxgBmFxkAEA3Beeo/sqISFBzrMLoxItrCFtiG2IvrHDUlil52eMp5gNAABuC0XYVxcvXyKEXNLWr5zkOESfryEz0okCQ/QAAG4NRdgnDzzwACMwC2NTznbUzUqynSajN/HbT2rM8zBEDwDg7lCEfbLph00ihrFyhqXpavsh+h+LNEzyPSQ0iGI2AADoDRTh3fvoo484nkwYeM/+urLF4xyG6L8+oTUsmEExGwAA9BIG6u/eb555RiGWxMiUgQkOQ/Q5ZVpTaDAZPohiNgAA6CWsCO/S2bNnzQyJ8w/bXVe6fGKPIXqtPgtD9AAAngErwrs0ceJEFSNJD46p8O+IsxuiP11tbNAxZNxoitkAAKD3sCK8GyaTSWsyysXS3JbzPa5E/3mexpSFIXoAAI+B5+u7kZiYKCfsgugEK2tIHWwbom/osBRfNgrTMUQPAOAxUIR341L1FUEgVZq6lY5Xov86T8PPSCdyGb1oAABwZ1CEd2zFihUsT7JiU8511M1MtJ0mozPxO0s1FgzRAwB4FBThHfvxpy0sw1g4w30ZDlei31LYSVJGkpBAitkAAOBOoQjvzIcffsgJ7MTIkQfqzy1Ocxii/7ZQa1wwnWI2AAC4CxifuDPPPvv/FCJJlFQRlKAO9LMN0R84qzEPCCXDMEQPAOBhsCK8AyUlJRaGHRwQvqfu1H0THIbo1+XrdBiiBwDwQFgR3oFp06apGPH4wKhKxyH60iuGRpMIQ/QAAJ7ICSvC8vLyJUuWxMbGRkVF9X1vbstgMOgsZrlYeqT53HVD9FrjwunEfpACAAA8hBOKUCwWL1my5K9//WtHR0ff9+a2EhMTFQK7IGo0IzGNiVN031/fbjl5xShMS6OYDQAA7poT3hodOnTo0KFDi4qK+r4rd3a5tkZORJWaulWOQ/Rf5Wm4mRMwRA8A4KFc+hmh1WptaWm5ePFi102ZTOYp76bed999Ip4sjEs52FDyamJ09/06E7/7lMb6xhTbpi3thOcpRAQAuFMMQ8KCaYegz6VFeP78+T179nzwwQddN+Vy+dGjR2UyD1hLbdu+VcSITFbdsgx/ici2HvzhRCcZG09CAv59+3Id+9I7fkEe8DcCANBfM/MvrSWjhtIOcmeMRqNWq+3lxnK5XCy+TdO5tAgTEhKefvrp1atXu/JB++7tt9/mBCYzKj6noWzDSttykOOFb09ojS/ca9vUwoUOUr7w+yAKKQEA7tA771+rtlhpp7hjcrlcpVLdfrtewxzh7b3wwvMKkSRCIpuVqA6wG6LPPqO1Rgwggz3j3V0AALghJxShyWTav39/QUEBx3H79+8/duxY3/fpPgoLCy0MOyQwYnfd6R5D9J8X6HSLMEQPAODZnPDWqE6n+/jjjwkhixYt+vjjj8PCwiZNmtT33bqJ6dOnqxhxWkBklX97bKhtiP7kZUOTSUTGjqKYDQAA+s4JRRgcHPz999/3fT9uSKfT6a2WMKnqcHPZs0sdhujX5WuNWTMwRA8A4OnwGeGtJCQkyHlm7sB4iczSY4j+VDWG6AEAvAGK8FZqGuoZhqnorF0xyWE5+GWuhps5kcikN/tFAADwFCjCm1q6dKmIJ4tix1RpGmck2E6T0Rr53ac01tmTKWYDAABnQRHe1I6dO1iG1Vk0901Q9xiiZ1JH2YboAQDAk6EIb+yNN97gBWZq1MhDDeVZqbYr0Vs5YUOh1pg1g2I2AABwIhThjf3hpZdkrChcIpuT5DhEf1ZriQwncRiiBwDwEijCGygsLLQwzLCgyJ01p+6fqLb/o88L9HpciR4AwIvgCvU3MD0zU8VKxgVEhgS2RQdLuu8vuWxoNrNkzEiK2QAAwLmwIuxJq9XqOatSKj/cfLbH1MS6PK0x614M0QMAeBMUYU8JCQkKQTQvIl4qsyTH2oboa65aTlUbhanjKGYDAACnQxH2VNvUQAg5r6lZMdlhOfhVfic/ZxKG6AEAvAyK0MG8efPEHLMoLqVK0zRjtG2IXmPg957WYogeAMD7oAgdZOdkM4RoTZ3LJ6jFDkP015hxo0mQP8VsAADQH1CENv/4xz94nkyLHnWosSJrXI8hep1xIYboAQC8EIrQ5uWX/yhnxQPE0nnJan+FbYh+/xmtNTqCxEZSzAYAAP0ERfhvubm5FsIOD47eVXvq/gkOQ/Sf5WkxRA8A4K0wUP9vs2bNVIkkY9XhwQGtUXZD9MWXDFc5CUm5h2I2AADoP1gREkLItWvXDJxVJVEcai5bfv0Q/eKZGKIHAPBWKEJCCElMTFQIorkD42WKnkP0Z+pMwpRUitkAAKBfoQgJIaS+pYkQUnbtyirHIfov8zq52ZOIVHKT3wMAAI+HIiSzZ88Wc8ziuLGXdc2Zo2xD9J0Gbt8ZDNEDALgXhUJx+43uBIqQHDp8kGFIp+lajyH6TSc0JC2RBKpv8bsAAOAyjKZTRPiMjAzn7tbXi/CVV17heJIZPepIY8WicbYvjrFywveFGtOC6RSzAQCAjUDImTOxMXFO37GvF+Ff/vJnGSsOEUnnpahVctvR2HtaY42JJLEDKWYDAIAunMCXNl1kWKakpMTpO/fpIjxy5IiFsCNDonfVlt6X4fAW6Lp8nX4RhugBANzCzprTBotJLZOpVKrbb32HfHqgfu7cuUpWkqwOD3Icoi+8qG8nMpI0gmI2AADo9vq5bC1vOZF9qD927rsrwvb2dgPH+cv8chrPrJjk8BJjXZ7OgCvRAwC4h5Kr1Rc7GqUMSU3tl6lu3y3CxMREP8LOjRip9OMSB8m776++ai6rxxA9AIC7eP18tpGz/u1vf++n/ftuETZebWEZ9sy1y6umOCwHv8jVcnMmE4lPv2kMAOAOBELaTbqDNWdFLHn22Wf76VF8tAinT58u5phFg1Jq9a2Zo2zfJtNp4LLPaqyzJ1HMBgAAhBCLVTh5RVfQVs0LZMH8Bf33QD5ahMeOH2MY0mHquH+CWsTaPgv8vqCTSU8iARiiBwCgbO9pzajAqO3VpbxI2LJlS/89kC8W4csvvywIzIzohGNNFxam2jrPYhW+L9TiSvQAAO5g0zH9YNVAhjBDBsX16wP5YhH+/R9/k7KiABE7f4x/jyF6PjaKRIdTzAYAAISQgio9scp21pcZGK4/hujt+VwRZmdnWwk7KjRmT+3pZRmOUxMYogcAcA8bj+omh8UbLCa1TO7n59evj+Vz50ZmZWUpWUmiakBQQGtkkG2IvqBK387ISOJwitkAAIAQUt1qrmgwt6mqtbyl+ODR/n4431oRtra2Gnk+SKE60HhmRc8r0euMi3AlegAA+jYc1S6ISblyrVnKkJSUlP5+ON8qwuTkZD/CzgwbEagWRsfYDdG3ms83mITJYylmAwAAQkiHjjtwVnPZoDHy1jfeeNMFj+hbRdjU1soy7OmOSz2Wg5/naqxzp2CIHgCAus0FnbNjRuU1VIhY8tRTT7ngEX2oCKdOnSrhmcWxKfXGq1PjbUXYoeMOnNFysyZSzAYAAIQQi1X4qUDLEBnPkMVZi13zoD5UhHkFeQwhV43tKyao7IfoN57oZCYmY4geAIC63aWahKDoHTWneEbYuHGjax7UV4rwpZdeEnhmZkxiblPl/LG2K9FbrMLGQq1pIa5EDwBA3w+5ukHKCCIIw4cMddmD+koRvvHGG1JWpGbZBWMdhuh3n9LwQ6JJFIboAQAoy6/Us1b5zvqzRlYoLi522eP6RBHu27fPyjCjBwzaVXvq/gkOb4F+XqDTZ2GIHgCAvu+P6SaFxRstZn+5Qi6X3/4XnMQnzpNcvHixkpUkKkODAlrCA2x/5fxK/TVWThKGUcwGAACEkEvN5soGc5vyipa3lB7OdeVDe/+KsKWlxcjzwQp1duOZ5T2H6LWGxRiiBwCg77tjmvkxyTWdLVKGSUxMdOVDe38RJiYmKoloVvjwYH9hdLRtrX2p2VzeaBEmjqGYDQAACCHtOu7QOd1lg8bIc2+99ZaLH937i7D1WgfLsifbL/ZYDq7P01jnYYgeAMDljCb7GTZCyOb8ztnRo/IbL4hZ4Ve/+pWL43h5EU6aNEnCk8WDkhuNbVNH2q410a7jDp3TcrgSPQCAi7W0c1V1Ywcruu+wWIWtJ7Q8kfICuW/Z/a5P5OVFeKKwkAik1di2YqKatfu7bizoJBNTiFp5818FAADnk+w8mJXqr5TZnpF3ndQkBcfsqi7lRcI333zj+kjeXIQvvPCCIAizBiUeb6qaP9bhSvSbirSm+Zn0ogEA+CSDiT144oHxtidkQSAbj2uj/cIJIfcMpXMhPG8uwrffelPCiFQss8jx1cfOUg0/LAZD9AAALsZk56YN87MfY8ur1EkFv221p42sUFRURCWV1xbhzp07LQybHDZ4d+3pZRkOrz4+z9NiiB4AwNV4Xr7r0Op0lf193x/TTRww0sJZA+R+MpmMSi6vPWfy/vvuU7GSeGVwYGBTj1cfnVI/Msp132IHAACEEJJ/KlpF7MfYLjabrzRZ2/0uawXL6aP5tHJ554qwoaHBKPDBCnV2w5n7J153JXoM0QMAuJzf1uxHMxyekDcc1cyNTqrTXJURZvTo0bSCeWcRpqSkKIl4Zvjw0EAyyvHVR2WLFUP0AACuVn5Z1tE+5R6HMbYj53UXDZ0GzvLBhx9SjOacIly3bt38+fPvu+++w4cPO2WHfdR67RrLssVtVSsm+dnf/3muxjJvKhGLaAUDAPBN8m3Zq9Mdxth+yNPMjU4oaKwSs+TRRx+lF80ZRfjNN9+8/PLLv/nNbxYuXJiVlXXhwoW+77Mv0tPTpTxZPCi52dQ+ueerDy03cwLFbAAAvqiljSmrzLIbYzNbha2FGrPAcgK/auXPKEYjTjlZ5t133/3zn/88Z84cQsiRI0c++uijN954o++7vWslpSViRtRsbFs1yeHVx4YCDZk0FkP0AAAuJtlxcNFYf4XUboytRJMSPGh3zWmBJevXr6eYjfR9RSgIwsmTJydM+Pcya8KECbQGQbr89re/5Tkya1BiQfNF+yF6k0X4obDTNG8axWwAAL7IYGQPFa5Kcxhj23RcG6kYwBASP+IeitG69HVF2NbWZjabg4ODu26GhIQ0NTXdbONz587l5OS8++67XTdVKtWWLVucOzjy/vvvyUQSPyJkpaodXn2c1JARsSQqzImPBQAAt8Xuy00frrQfY8u9oPNjVNvrzhhZIScnR6vV9t+jy+Vysfg2TdfXIlSpVAzDGAyGrps6nU6tVt9s48GDB8+fP3/hwoVdN2UyWUhISB8D2Pvpp5+shB03IG537el1iyO77xcEsr5Aq3/yPic+FgAA3B7Py3YffnhZgP193x/TjQ9N2nilONBP6dwWuDt9LUKZTBYWFnbx4sW4uDhCyKVLlwYNGnSzjeVyeVxcXGpqah8f9GZWrVqlFInvUQYHBjX1ePWhlSsxRA8A4Gq5pTH+jP0YW1WTqbrZ2u53SctbyvNKKEbr5oSzRletWvXxxx8TQrRa7TfffLNq1aq+7/Mu1NTUmAQ+xC9gb/2p64botYZF+E41AABX89ue/VjPIXrt3Kjkek2bjGGGD6fzLds9OKEIX3rppYqKilGjRo0YMSItLW3JkiV93+ddGDt2rJIRzxwwLDKYjY9yePVR1coJE1KopAIA8F3nL/l1Xpt0j60I27Tc0XLdBcM1A2f5178+ohjNnhPGJ8LCwoqKiioqKlQqVVRUVN93eHfaNBo1Ky1qr/zZTMflYK7WMn8ahugBAFxMvi374fFq+2vRb8rrnBed8GP1GTFDVq9eTS+aA+d8swzDMPfccw/FFkxLS5PyJCsmqd3S89XHsXItd28GrWAAAD6quY05V7VgjO1bTYwWYVuhxkRYThAeeOBBitF68JLvGi09VcowTLPx6vKJDq8+NhR0CpNTMUQPAOBiku05S1N7DNF3jgmJ3V19mrBk3bp1FLP14A1F+MwzzwgCM3tQYmHLpXmOrz42F2nMCzLpRQMA8EkGI3ukaIXjEP3mXP1AxQCGkNHxoyhGu543FOG//vU/EoaVE37ROHWPVx/knsEkIpRiNgAAH8TuPT5hhDLM33YayrEKnZJRbq87Y2T5EydOUMx2PY8vws2bN1sJOyZ8yJ7aM0vSHV59fHFCp8+6l2I2AABfxPHSXYceyXC4Ev3Go7q0kBFWzhqk9L/tV724mHuluQsPPvigUiQe4RcUeN2rD61CReKHUMwGAOCLck/GBYlGDLR9fWZlk6mm1drud0krWCvycilGuyHPXhFWV1ebBD5MGbS3/tSKSQ6vPtbl6fSLZ9EKBgDgs/y2X3cl+iPauVHJDZo2GcMMGzaMVrCb8ewi7BqizwwdEhXi8OrjUov5UpuVZM1qqpQAACAASURBVCRRzAYA4IvOX/LTdE4cYSvCVo31WLm2wtBu4CyffPIJxWg349lF2KHViFlRUXtljyvRt3ZaJZGhRIQhegAAl1L8tH91uspxiF4zLyapuPGimCEPPuhG44PdPLgIx44dKxXYrJikDmvnhBGYFAQAoK25jZRfmp9iO2/RaBF2FGkMAsMR4aGHHqYY7RY8uAjPlJ0hhDToW1ZMUtq/+gAAACqk2w4scxxj21HcmRo6eG/NacIK7vm+KPHcInzyyScFnpkbm1TcemVu8k2vgAgAAC6iM5AjRcvHOQ7R5+nC5SGEkKQE9z1pw1OL8NPPPpUwrETgFqc5vPoAAAAq2H3HJ49U2Y+xHS3XqVn19rrTZlbIy8ujmO3WPLJCNm3aZBXI2PAhe2rPLEvHchAAgDaOk+4+vNpxiP77o7pxwcOtHBesdrshenvum+wWHvr5z5UiyTC/wKB7lCFqj/wrAAB4leMnhwSLhkfYxtjKG0z1bVybrEorWE8VFFCMdluetyKsrq42ESFcFbS37tRyxyF6AACgwm/7gR5D9N8d1cyJSm7WdsgZJi4ujlKuXvG8IhwzZoySEU8LGTJogMOrDwAAoONslVLbmTHcVoQtnda8C/rzuqsG3vL5+vUUo/WGhxUhx3Edep2YFRe1XVgxCbODAAD0KbZlP5Jx/ZXoE0ubL4sZYcWKFfSi9YqHFeHYsWNlPLMoJukap0kfjiIEAKCtoZWUX5qb7HAt2J3FWp1ALAL/i188QTFaL3lYEZaVnxMIqdM1r5yswhA9AAB10u0596c5XIl+e3Fn2oAh+2rOMKzwwQcfUMzWS55UhGvXriU8M29QUvHVy3MwRA8AQJ3OwBwvtr8SPS+QTce1obJgIggpSSkUo/WeJxXh55+vkzCsWLAuHe8vl2A9CABAGbv36OSRKvsxtqPntUHigK01pWYRyc11u0sP3pDHFOGGDRusAkmLGLa37uzSdH/acQAAfB7HSfccXZ3uOER/TJcaPEwQhNCAQJGHXALIY6bRH3lktZ9IOliuDhypDFF5xsEFAPBmR0uGhYiG2Y2xna83NbULB7RVWmI95cbfqdaDZ6wIKysrzQKJUAXurT99/0QM0QMA0Oe3PXvNdUP0syOTWvUeMERvzzOKcPz48SpWvCQycYA/OwxD9AAA1J2pVBq06cNsRdjcac2v1J/TXtVx1m+++ZZitDvlAUXIcdw1g07MilNDY1VyvCkKAECfYlv2mgkOQ/Q/5HbOj0kqbbksYYQlS5bQi3bHPKAIk5OTZTz7zKgZYsYD0gIAeL+GFrbyyrwk29SEwczvKNZe43iLwP/yl09SjHYXPKBayisrGEIeGTaRdhAAACCEEOm2nPvGqWV2Y2zbizQZYUOya88yrPDuu+9SzHYX3L0IH330UYZnVg6fECBV0M4CAACEaPVMbsnyNNsYGy+QTbnaIGkQQ0jqmFSK0e6OuxfhV19/KWbYp0ZMpR0EAAAIIUS05+jUeJX9GNvhc9owWdC22lNmVjh69CjFbHfHrYtw3bp1VsIOC4kerA6lnQUAAAixcuLdRx66bog+JWgoL/ADAoM8ZYjenlsP1D/55C+Z4CCFTE47CAAAEEIIOVo8Ilw6LNxuiL7O2NpB9usu6ASuoriYYrS75r4rQovFwnMWfvhI2kEAAODf/LZn97gS/Yaj2lmRiVf11+QsO3DgQFrB+sJ9V4RarVagnQEAAGxOV/ib9eOHBnbf0dxpLajSDw9u1XHWn37YRDFaX7jvihAAANyKYlv2Ixkqxm6IfuPxzoWDUs60XpEwQlZWFr1ofYIiBACAXqhvZqtq5toN0etN/M4SbbvVahH4p556mmK0PkIRAgDA7Um3HVie5i8VOwzRTwgfeqD2LMOSN998k2K2PkIRAgDA7Wj1TF7p/ePsrkTPkx/ytIGSAEEQ0tPGU4zWdyhCAAC4DdGuI5mj1MF2Q/SHzmnD5cHba09bROTIkSMUs/UdihAAAG7JYpXsPfpwRs8h+qSAOF7gw4KCaeVyFvcdnwAAALdwtHjEQOngAdLuO87VGds6Sba+SidwFUVFFKM5BVaEAABwK37bDzya7jBE/+0R7cyBiVf11xQiTx2it4cVIQAA3FxpeYBFnzbENkTf2GEpvKhvD27R8dZtmzZTjOYsWBECAMBNdV2J3mGIPlezICa5rLVaSoQFCxbQi+Y0KEIAALiJ+mb2cu3sRNtpMnoTv/ukto2zWAT+6Wd+QzGaE6EIAQDgxmQ/Za9MU9sP0W8t7JwUPiyntoxlyeuvv04xmxOhCAEA4EY6taTg1P32V6LnyeZ8nb80gAhk4oSJFKM5F4oQAABuQLTn6IzR6kA/2xB9Tpk2Uh6yrfqURURycnIoZnMuFCEAAFzHYpXsPXb9legTAuIEgQ8PDqGVqz9gfAIAAHpijhTFD5TG2Q3Rn642dmjIXl25jnBVJSUUszkdVoQAANCTYvuBNY5Xov/+mHbWwKQOo1YhEoWFhdEK1h+wIgQAAEel5QGcIXVwUPcdDR2W4suGjqAmHW/d9eNWitH6A1aEAADgQLF1/2M9huiPaxbGJJ+7Wishwpw5c+hF6xcoQgAAsFPTKLpSNzPBdpqMzsTvPqm9arWYee73v/8PitH6iXOK0GQyFRYWetPZtAAAvkm27cDK8Q5Xov/pROeUgcMP1J5lGfLKK69QzNZPnFCER44c8ff3X7Ro0cKFC/u+NwAAoOaalpw4fZ/dleg5XvghT6sSq4nATJ48mWK0/uOEIkxJSWlqatq2bVvfdwUAABSJ9hyZmeA4RH9WG6scsK2m1CISsrOzKWbrP04oQrVaHRgYePvtAADAnVmskn3HHxzvMES/8Zh+VEAsEcjA0AG0cvU3l45PtLS0bN++va6uruumWq1+4oknWPbGZWyxWFwYzWn0HeYDeztppwAAuL3OFpP9TeZw4ahIhyH60isGrY7doy/XM3x5fr4nPi2LRKKbtUy3XhVhQUHBmjVrrr9/69atQ4cO7X0gs9nM83xbW1vXTa1WazabJRLJDTfmOK73e3YXkQMMUyfs5gXaOQAAbk+YQEhs5P/dEBTbsh+d0+M71bTTIxK+v1KiELHBwcGe+LR82xYkvSzCxMTELVu2XH9/TEzMHQWKioqaPn366tWre7OxXC6/o527BaWCeyCLdggAgDt3sjxQMI+Js32JaH27pfSKsSOwUcdbDuza45HPyb3TqyJUKBQjRozo7ygAAECLYtv+xyaqrhuiT9l4+aSECJmZmdSS9T8nfEbY0dHx/PPPt7a2ms3mtWvXhoaGvvrqq33fLQAAuEhNo7imYeZ90d136Ez8nlJt+kCTmef+849/oBjNBZxQhBKJJDU1lRDS9b07arX6dr8BAABuRLYte9V4tURkWw/+WNA5LXLEntoyEUv+67/+i2I2F3BCESqVyieeeKLv+wEAAAquacmJM8ueti0HOV7YnK/NDFcRQqZNy6QWzFXwXaMAAD5NvOvwrER1gN0QffYZbZwybFtNqYUle/fupZjNNVCEAAA+zGwR7TvWY4h+03F9vH8MISQqLJxSLJfC9QgBAHwXc+hEQowiNtQ2RH/yskGvF+3Wl+sJV1NaSjGby2BFCADgqwRBsT3n0Z5XotdNDx/dadT7iSU+8vWZWBECAPiqkvMhrGVMnKL7jvp2y6lqQ3tAvY43H9yxn2I0V8KKEADARym27X90gsNy8LtjmqxBYy60N0gIM3XqVFrBXAxFCADgk6obxLWNM0bbTpPRGPi9pzQNJoOJt/7xj3+kGM3FUIQAAL5ItjX7gXTHIfoTndMjRx6uPy9iyX/+539SzOZiKEIAAN9zTUOKzi4Za/siMCsnbMnXykVKQoQZ0++lGM31UIQAAD5HvPPwnKSeQ/RD1eHba0qtrLBr1y6K2VwPRQgA4GNMZtH+4z9Pd/he6E3H9SPU0YJAoiMib/Z73grjEwAAvoU5dCIpVhEdbLsoevElg8Eg2q0vNzDcyZMnKWajAitCAABfIgiKHTmPpvcYotdmho/WmHRKiTQgIIBWNFqwIgQA8CXFZaEia3KsbYi+5qrldI2x3b9Ox1uP7fGVIXp7WBECAPgQv63ZjzkO0W883pk1KKWyvUHKkAkTJtAKRhGKEADAZ1xpEDc0TR/lMES/77S2wWQwC9x//def6CWjCUUIAOAr5Fv3PZCuFtsN0W8puDY9cuSR+vMilrz44osUs1GEIgQA8A3tnaT43LJU/+47rJzwY4FOKvITBDJr5myK0ehCEQIA+ATxrsNzk/1VctvT/r7TmmHqiB3VpziRsH37dorZ6EIRAgD4AJNZlJ3b40r0G4/ph6mjCBFio2Jo5XIHKEIAAO/HHCxIiVVE2Q3RF13Um02SXQ3nDAzvg0P09lCEAADeThAUOw+uyehx6UFdZvgondmgkkhVKtXNftUXYKAeAMDbFZ4dILYmDbIN0VdfNZ+vM7Wr63S8NW9/DsVo7gArQgAAL+e3recQ/ffHNAtiUi51NEqIMH78eFrB3ASKEADAq12qFTc2Z8bb3vzsNHDZZ7S1Jp2J5/72t79TjOYmUIQAAN5MtjX7oR5D9Pmae6Pij9WXi1jh2WefpZjNTaAIAQC8V9s15uT5xXZD9BZO2JyvEbEK3reH6O2hCAEAvJZk56H5KWr7Ifq9pzQjAyN3XCnlRMKGDRsoZnMfKEIAAC9lMrM5+Q+M73kl+qGqSIZhhsTE0srlblCEAADeiTmQPzZOERlkG6I/UaUXzNKd9WUGhispKaGYza2gCAEAvJEgyHfkrMlwmJT//ph2Sli8wWpSy+RKpfJmv+prMFAPAOCNTpwJl/EJMfLuO6pbzeX15g5VrZazFOUcoRjN3WBFCADghfy2Zf9iosNy8LuuIfprTVKGjBkzhlYwN4QiBADwOhdrpS2t0+Jtb352GrgDZ7XVRq2Rt77++hsUo7khFCEAgLeRb93/ULpKxNqG6Dfnd86KGpXbUCFiyNNPP00xmxtCEQIAeJe2a6S0PGuM3RC9VfixQMuwcp6QhQsWUozmnlCEAABeRbLj4MIxjkP0pzXxAVE7qkt5Vti8eTPFbO4JRQgA4EUMJvZA/qo0xyH6Y/rBqoEMIUNjB9PK5c5QhAAA3oPJyRs31M9+iL6gSk+ssh31Z3El+ptBEQIAeAuel+849Ei6w9TExqO6yWHxRqtZLZPL5fKb/aovw0A9AIC3OHEmSimMdhyir2gwX1VWaznLycPHKUZzZ1gRAgB4CcXW7EczHL44bcNR7YKYlOrOZinDJCUl0Qrm5lCEAABeoapGfvXq1JG290U7dNyBs5pLBo2Rt7777rsUo7k5FCEAgDeQ/7T/oQwVa/ekvrmgc3bMqPyGCjHLrF27ll40d4ciBADwfC3t5FRFjyH6nwq0hMh4hrz6yqsUo7k/FCEAgMeT7DqUleqvlNme0neXakYFRu24clIg/O9+9zuK2dwfihAAwMMZTOzBgh5Xov8hV+cnUfMMeevNt2jl8hQoQgAAz8YcyE0b6hceYBuHy6/UCxZZTt05gXBPPfUUxWweAUUIAODJeF6+89Dq9B5XoteFK0IEhhQUFNLK5UEwUA8A4MnyT0eryOho2xD9pWZzRb1JZ7k8KnF0cnIyxWieAitCAAAP5rdtf48h+u+OaeQimcCS4uJiWqk8C1aEAAAeq6pG1t4+5Z7o7jvq2y05ZRpOYI/n5lLM5VmwIgQA8FTyH/etTlfbD9Gvy7lm5YR7Ro8aO3YsvVweBkUIAOCZWtqZsxeyxtqmJmrbLAfOdnIsizdF74jTitBoNFqtVmftDQAAbk2y4+Cisf4Kqe1p/OP97YRhCwtxpuidcUIRbty4MT4+PigoKDAwcMmSJVevXu37PgEA4FYMJvbQCfsr0ZdeMRwt18QnJCYmJlLM5YmcUIQymeyLL77Q6XSNjY16vf6FF17o+z4BAOAW2P256cOV3UP0DR2W579psAqioqIiusE8kROKcNGiRWlpaSzLqlSqpUuXlpWV9X2fAABwUzwv23Xo4fH/nprQm/jnvmi0WNmysrN0c3koZ45PCIKwZcuWzMzMm23AcVxLS8vFixe7biqVyvDwcCcGAADwCXmnYvyZUdFyQggvkD9939LcaR07Pn348OG0k3mkXhVha2vrhx9+eP39Dz/8cFxcXPfNV199ta6ubsuWLTfbT1lZ2e7duz/44IOumzKZ7OjRozKZ7IYbazSa3mQDAPA1ftv2P/Z/Q/Tv72o7V2/mGNGuXbvu9GlTp9PxPM8wTD9kdBdyuVwikdx6G6etCP/5z3+uX7/+0KFDSqXyZtskJiY+88wzq1ev7s0OcQ4qAMANlF/y67w26Z4oQsiOIk32aYPBItTV16nV6tv+ag8MwyiVSu8uwt7oVRGGhoa+/PLLt9jgk08+eeuttw4ePBgZGemkYAAAcAPybQceHq9mGVJ6xfDh3naDmf/Fk2tDQkJo5/JgTlgRfvXVV0899dR7771XWVlZWVmpUCgmTZrU990CAEBPLW1MWeWC/xdT327547ctejOnCgp47733aMfybE4owqampilTpnz//fddNyMiIlCEAAD9QbI9Z2mqPy+Q//iyxcwxVhHT2tJKO5THc0IRPvfcc88991zf9wMAALdiMLKHi+5fO/AP37TojKyJWC9eqKSdyRvg6hMAAJ6B3Xd8wgjlN0c6a1sFjdnyv//78aBBg2iH8gb40m0AAE/A8dKdhyIVwuEyY4fJsmjZkocffph2Ji+BFSEAgCfIKx3gx+ws1pssQsSgqG+//ZZ2IO+BIgQAcGM8T85dlOaWkGMl7RbOyjMimaT7+7nAKVCEAADuh+dJxRXp8SImtzRQzswYLssRMVqL2MLyBp2OdjhvgyIEAHAbPE8qrsiOFwnHTgb5sfNHyec8MiAyUPLb9c0CLzEz5itVWAs6n9cWocDxRGegnQIAoBd4npy/JD9exJecjwn3WxAvmf5ERPclll776Wp9G7lmtnz//XdRUVF0k3ol7yxCtUJkrW6Q/fq/aQcBAOiVwZGKhSNl056KDlGJ7O//5ui1vHJzp8n68JqHFy9eTCued/POIhwZKTvwQhztFAAAfZJbof/mSKfBwscNH/rRRx/RjuO1MEcIAOCOLrWYX/mhVWfmRHIpLnjer7xzRQgA4NHatNxz65vMVsKJGCMuztrPsCIEAHAvZqvw4tctVqvYzAjNjY2043g/FCEAgBsRBPKPLVdbO4nOat2xfXtQUBDtRN4Pb40CALiRzw92lFyyaMzc2l/9ctasWbTj+AQUIQCAuzhYpvshT2Ow8CMTR7/99tu04/gKvDUKAOAWyhtMr/90VWfmpEpFcXEx7Tg+BEUIAEDfVY31xa+aTVaBEzMdHR204/gWFCEAAGUmi/DC1y0Wq8hMhI6rbbTj+BwUIQAATYJA/rb5aruG0XHWI4cPK5VK2ol8DooQAICm/81uP1Nt7TBann/xhYyMDNpxfBHOGgUAoIPjhZ8KO7cX6nQWLjV9/J/+9CfaiXwUihAAwKUsVqHwov7gGeORcm24wl9rsqoC/Y8cOUI7l+9CEQIAuILZKhRU6Q+eNh6t0N4TOGBe1ORHJkU/fOgzXsK2tLbSTufTUIQAAP3IZBFOXNQfPG08Wq65Jyh86aDMNxaOjfALaDPppuz4h4nhDRod7Yy+zsOKsElj+qnwGu0UAAC3Z+WEoipL4UXduPCYJdFT31mcHCpXdf2RhedWHvy43WIoLCgQiz3sedj7eNJ/gPigyMkhifXlPO0gAAC3xzLs8pDhn6UkBkr9evzRbwq+O9/e8Oe//Dk5OZlKNrDnSUUY6Rf49vif004BAHD3rpq0/yw7sPVS8cRpU/7jP/6DdhwgxLOKEADAQzUbOrfVlH5bd+ps8xUxIf4hQfv27aMdCv4NRQgA0F/q9R17a89ubCwrbrggGjbcEBdOrtWzVktDQwPtaGCDIgQAcLJaXfu26tKv605WdTSyoxKMg4IkzABj5XnCMDKB0xsttAOCAw8owpr2xj8Xb6WdAgDg9iwCt6+lqlbXRpKSjQlDJZVW0+kShmFY3iplJcuWLfvyyy9pZ4Se3L4Ig4MbFs57l3YKAIBe0WqZmHjJhfPW4gKGEEbgFBL5Y4899u67eBpzX+5ehOKD+8QBKtopAABuj9MZeb2R8JyYEfwUqj/84Q+/+93vaIeC23PfIlQoFGKRiDHqiVFPOwsAwO2JCatQBvz1r6+uXbuWdha4A+5bhHK53OCBHylrNBq1Wk07hWfAseo9HKvew7GCO4XrEQIAgE9DEQIAgE9DEQIAgE9DETpTa2vrp59+SjuFZzCbze+88w7tFB7j9ddfFwSBdgrP8D//8z8ajYZ2Cs/w9ddf19bW0k5BH4rQmaqqqr7++mvaKTxDW1vbBx98QDuFx3jttdfMZjPtFJ7hk08+qa+vp53CM2zcuLGsrIx2CvpQhAAA4NNQhAAA4NNQhAAA4NMYV34Cn5GRcenSJZXKa78yzWw2t7a2RkZG0g7iAXier62tHTRoEO0gnuHy5cuxsbEMw9AO4gFqa2sjIiLEYvf9thD30djYGBgYKJfLaQfpRw888MBf/vKXW2/j0iLs6OhobW1lWW9ehppMJplMRjuFZ8Cx6j0cq97Dseo9XzhWAwcOVCgUt97GpUUIAADgbrx5cQYAAHBbKEIAAPBpKEIAAPBpKEIAAPBpOMPYOS5evFhYWNjR0fGzn/3sZtdC27Jly5EjR2JiYh5//HEvniG5LUEQNmzYcOLEicGDBz/++OPXn7p97Nixs2fPdt/8xS9+4d1nGveQn5+/adMmtVq9Zs2amJiY6zdoaGhYt27d1atXlyxZMmXKFNcndB/V1dXr1q3TarXLly8fP358jz81m82ff/55982UlJTrt/EROp2upKSkvLx86NChmZmZN9ymvr7+s88+6+joWLJkyeTJk10bkDIfen7pP3V1dePGjfvoo4/Wrl3b2tp6w21ef/313//+98OGDTt06NCsWbN8+WTdP/7xj6+++urw4cN37NixZMmS6zf47rvvvvzyy4v/x/UJKcrOzp4zZ05ERERbW1taWlpLS0uPDTo7O9PT069cuRITE7N06dKtW7dSyekOmpqa0tLS2tvbIyIiZs+enZOT02MDg8Gwdu3aqqqqrn9IbW1tVHK6g+eee+6Xv/zla6+9Zv/KwF5HR8f48ePr6uqio6MXL168c+dO1wakTYA+4zhOEASj0UgIuXjx4vUbmEymsLCww4cPC4JgsViio6P379/v6pTuobOzU61Wnzp1ShAEg8EQFBRUWFjYY5unn376T3/6E4109M2aNevNN9/s+jkrK+uvf/1rjw3ef//9qVOndv38ySefZGRkuDSfO3nllVcWLVrU9fMbb7wxZ86cHht0dHQQQqxWq8ujuZ2u56gXX3xx9erVN9zgnXfemTFjRtfP//rXvyZPnuyybO4AK0InuO0bd+fOndNqtZMmTSKEiMXi6dOnHzp0yCXR3E5JSYlSqUxMTCSEyOXyyZMn3/BQFBcXv/baa999913XywsfIQjCkSNHZs+e3XVz1qxZ1x+cQ4cO2W+Qn59vMplcmtJtHD58eNasWV0/z5o16/Dhwzfc7P3333/vvfdKS0tdGM3t3PY5qse/q+PHj1sslv7P5S5QhK7Q2NgYGhra/W8xPDzcZy8T09jYOGDAgO6bNzwUUVFRkZGRHR0db7311pgxY65du+bajNS0t7cbjcbu4xMWFtbQ0NBjm4aGBvsNBEG4fhsf0eNQGAyG9vZ2+w0Yhpk9e3Zzc/OZM2emTJny3nvv0YjpGXocTJ7nm5qa6EZyJZws0yutra0RERHX379hw4b777//tr8uFoutVmv3TYvF4t1fa5SUlHT9Rc5+9atf/fOf/xSLxRzHdd95w0Px/PPPd/3A8/yECRM+/PDDF198sV8DuwmJREII6f6nYrVapVLp9dvYb0AIuX4bH2H/v1XXD10HsJu/v/+ePXu6fl62bNnSpUvXrl3rs4fr1nz83xWKsFdCQ0Ptm+xORUZGtra2dn+tX11dXXJysvPSuZ1Tp07d7I8iIyMbGhp4nu9aH3edZ3SzjVmWnTBhgu+cL6NWq9VqdV1dXVRUFCGkrq5u4MCBPbaJiorqXkPX1taKxeKwsDBXB3UP9oeirq4uICDgFidjT5w40WAwNDQ0xMbGuiqgJ+lxMKVSaWhoKN1IroS3RvtReXn5+fPnCSEjR46Mi4v76aefCCHt7e3Z2dlZWVm009Exbtw4hUKxf/9+QkhDQ8Px48cXLFhACGlqasrNze3axmAwdP+wf//+hIQEWmldLysra+PGjYQQjuM2b968aNEiQojFYjlw4EDXYcnKyvrxxx+7Pr/ZtGnTvHnzfPYyC1lZWZs3b+Z5nhCycePG7v+n8vPzu94utv+Aedu2bQEBAdHR0VSiuiez2XzgwIGuz5izsrK2bNnS9e9q48aNCxYs8KmZJZw16hyZmZljx44lhCQkJKSmphoMBkEQHn/88UceeaRrg82bN4eGhj7yyCPx8fFr1qyhGpay9evXh4WFrVmzZujQoc8++2zXnV999dXgwYO7fo6JiVmwYMGDDz4YExMzY8aMroPpI8rKysLCwlauXDllypTx48fr9XpBELo+ramsrBQEwWw2T506NT09/cEHHwwNDS0qKqIdmRqdTjdu3LipU6euXLkyPDz83LlzXffHx8d/+umngiC89957iYmJP//5z2fPnu3v7//dd99RzUvTN998k5qaGhERERISkpqa+u677wqCUFdXRwi5fPmyIAhGo3HSpEkTJkx44IEHQkNDT548STuyS+HqE85x8uRJ+4++xowZw7Js17+wwYMHd9156dKlvLy8mJgYXxtWvV5FRUVRUdGQIUPS09O77rl69WpdXV1SUhIhpKam/P2wTgAAAOZJREFUpri42Gg0Dhs2LDU1lWpSCtra2vbv3+/v7z9jxoyuz2msVmtxcXFycnLXW+tWqzUnJ6e9vT0zM9Nn3xftYjabs7OztVrtzJkzg4KCuu48ffr0wIEDQ0NDzWZzYWHhlStXAgMD09LSfOq9vh6am5tramq6b0ZERERFRVkslpKSkpSUlK5/ZhaLJScnp6OjY/r06fZntPkCFCEAAPg0X3oXGAAA4DooQgAA8GkoQgAA8GkoQgAA8GkoQgAA8GkoQgAA8GkoQgAA8GkoQgAA8GkoQgAA8GkoQgAA8GkoQgAA8Gn/H92qpQh3+0tPAAAAAElFTkSuQmCC",
"text/html": [
"\n",
"\n"
],
"image/svg+xml": [
"\n",
"\n"
]
},
"metadata": {},
"execution_count": 7
}
],
"cell_type": "code",
"source": [
"plot!()"
],
"metadata": {},
"execution_count": 7
},
{
"cell_type": "markdown",
"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$.\n",
"Note that `A^s \\ W` triggers the computation of the H-representation of `W` and `A_W` is H-represented."
],
"metadata": {}
},
{
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": "αo (generic function with 1 method)"
},
"metadata": {},
"execution_count": 8
}
],
"cell_type": "code",
"source": [
"function αo(s)\n",
" A_W = A^s \\ W\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"
],
"metadata": {},
"execution_count": 8
},
{
"cell_type": "markdown",
"source": [
"We obtain $\\alpha \\approx 1.9 \\cdot 10^{-5}$ like in [RKKM05]."
],
"metadata": {}
},
{
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": "1.9190700000000013e-5"
},
"metadata": {},
"execution_count": 9
}
],
"cell_type": "code",
"source": [
"α = αo(10)"
],
"metadata": {},
"execution_count": 9
},
{
"cell_type": "markdown",
"source": [
"The scaled set is is the following:"
],
"metadata": {}
},
{
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": "Plot{Plots.GRBackend() n=1}",
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO3dfWCN9f/H8etsZ3e2uduaMVMi0jdSM/ezkajY3EVCjciSm0gligqlkNnmviG35SaTUSrMjORmyE0iCV8bY3fs9uzsnOv3x+mnfZndnp3Pdc71fPy1c7pc16urdV7e57o+52hkWZYAAFArO9EBAAAQiSIEAKgaRQgAUDWKEACgahQhAEDVKEIAgKpRhAAAVaMIAQCqRhECAFSNIgQAqJpFi3DNmjXHjx+35BEtT6/Xi45gNThXZce5KjvOVdkVFhaKjqAIFi3C3bt3nzx50pJHtLz8/HzREawG56rsOFdlx7kqO86VCW+NAgBUjSIEAKgaRQgAUDWKEACgahQhAEDVKEIAgKpRhAAAVaMIAQCqphUdACJlZ2d///33x48fl2XZ8kcvKChwdHS0/HGtEeeq7DhXZWeN58rV1XXKlCn29vZm3CdFqEa3bt2KjY1dveHbhL17XJq2z2zQXrYz528VAFQRh+9nhYWFeXl5mXGfFKGKZGRkxMbGrli/6df9+5ybBdxq1kOavjDf/QHRuQCgrLR7Isy/T7PvEUqTlpa2Y8cOU/85PBqQ/UQ/adZXOpfqonMBgCJQhLbs4MGDY96ecubUCW2LZ3NavCyFrNU5VhMdCgCUhSK0WX/99ddzIX1vhcyShmzTOTiLjgMACkUR2qbbt28/07N3zrPvS+0Gi84CAIrGOkIbZDAYeg0YlPxg58LAMNFZAEDpKEIb9PrY8Ueu63V9PxcdBACsAG+N2ppl0cu/3r4r5+19kh3/cQGgdLxW2pT4+PgJ703NfWuPVK2m6CwAYB0oQttx7ty5nn0H5L66TvJqJDoLAFgNrhHaiPT09KefC84N/lhqGig6CwBYE4rQFuj1+h59B6Q+1tvYcbjoLABgZShCWzBy9LiTtx11IdNFBwEA68M1Qqs3d174ph/35b69T+IbJACg/ChC67Zz586PPvsi5+19Eh+iDQAVQhFasbNnz/YfMjRn5Gaptq/oLABgrbhGaK3S0tK6Ph+S22e29HAb0VkAwIpRhFZJr9d3D+mX+sRAY5uXRGcBAOtGEVqlYSPfOFvoUdBjquggAGD1uEZofWZ8MmtrQmLuhD2SRiM6CwBYPYrQymyJifls/sLcd/dLTq6iswCALaAIrcmJEydefnVk7hvfSbV8RGcBABvBNUKrcf369Wd69MobECE19BedBQBsB0VoHfLz87sF97nVboTs3190FgCwKRShFZBleVDo8AtaX333SaKzAICt4RqhFZj64cc/n/w7782fuE0UAMyOIlS6TZs3z/9yVc47CZLWSXQWALBBFKGiJSYmDgsbkzP2e6l6HdFZAMA2cY1QuZKTk7uH9M0ZtFiq31x0FgCwWRShQuXl5XXt0ftWwBjpiZ6iswCALaMIlUiW5QFDhl6q3qyw63jRWQDAxnGNUInefm9K3NnkvHE/iA4CALaPIlScVatWL1m9Iffd/dwmCgAWYIYilGX5yJEjJ0+edHZ27tq1q7e3d+X3qVq//vrrGxPeyZ2wS3J/QHQWAFAFM1wjnDx5cmho6OHDh2NjY5s2bRofH1/5fapTWlpaz74Dcocsleo+KjoLAKiFGSbCcePGzZo1S6PRSJI0efLkWbNmBQYGVn63KvRS6PDsJwdILXqIDgIAKmKGibBevXqa///or9q1axuNxsrvU4UWLlr8yx9XdMEfiw4CAOpizptlUlNT58+fHxUVdb8NUlJStm/fnpSUZHro5OQ0evRoe3t7M2YQTq/X6/X68v6ps2fPvvvBh7kT90pax6pIBQA2o1wvs/b29nZ2pYx8ZivC7OzsXr16vfDCC3379r3fNgaDIScnJyMjw/TQ2dnZYDBobOuDpI1GY3ln4oKCgn6DXskPmSHVeaSKUgGArZBlWS77y2ypLSiZqwhzc3ODg4Mfe+yx+fPnl7BZvXr1OnfuHBoaapaDKlNBQYGTU/mWPUybPvOay4PGjq9WUSQAsCEaR0fH8r7MlswMRajT6fr37//ggw8uXbrUxsY7Czh+/PiCpdG5kw+LDgIAKmWGm2WmTJmyZ88erVY7atSosLCw999/v/L7VAmdTtdvUGhev7lSDRZfAoAYZpgIX3rppbZt29556OrqWvl9qsQHH358w62h7D9AdBAAUC8zFGGrVq1atWpV+f2ozbFjxxZ9uTL3/SOigwCAqvHtE2LodLoXBofm9Q/nG3cBQCyKUIz3PpiWUrOZ3Kqf6CAAoHZ8+4QAhw4dWrpidd4HR0UHAQAwEVqcTqcb+Mrw/AHzJXcv0VkAABShxb075YMbns1lv/t+/g4AwJJ4a9SiDh48GL16fe77vCkKAErBRGhRI8ZMyO07W3LzFB0EAPAPitByYmNjL6dlS61eEB0EAPAvitBCZFmeNG1GTs/pkoZzDgAKwouyhcTExFzNNkhP9BQdBADwP7hZxhJkWX77/Y+ynp8p8e0cAKAwTISWsHPnzlS9vdT8OdFBAAB3owgtYfrn87KDxotOAQAoBkVY5U6fPn3yzFmZm0UBQJEowio34/MvdIGjJK2j6CAAgGJQhFXrxo0bsdu2GToOFx0EAFA8irBqzYuIktsMlNw8RAcBABSPIqxCeXl5i5Ysy+/0huggAID7ogir0FerVkmN20neTUQHAQDcF0VYVWRZ/vSLyKzAcaKDAABKQhFWlR07dtySXKQmAaKDAABKQhFWlRmzw7NYRA8AikcRVolTp06dPntOatVPdBAAQCkowiox/bO5uqDRkr2D6CAAgFJQhOaXnJy8PTbW0GGY6CAAgNJRhOY3P2qh1HYQi+gBwCpQhGaWl5e3ZFl0ftAY0UEAAGVCEZrZmrVr5cbtJa9GooMAAMqEIjQnWZbnLVyaHcQiegCwGhShOcXGxt7WuEmPdBQdBABQVhShOU2fHZ7dZYLoFACAcqAIzebkyZNnz1+Q/PqKDgIAKAeK0Gw+njWHRfQAYHUoQvNITk7+/vvv+SZ6ALA6FKF5zIuIktoNlqrVFB0EAFA+FKEZ5ObmLv1yeX7gaNFBAADlRhGawfLlK+QmAdIDD4sOAgAoN4qwsoxG4+fzF+SwiB4ArBNFWFnbtm27ZecmNW4vOggAoCIowsqaMTs8u8tbolMAACqIIqyUxMTEcxcvS359RAcBAFQQRVgpM2fPyw98Q7LTig4CAKggirDikpKSdu7caej4quggAICKowgrbu78SLndy5JLDdFBAAAVRxFWUE5OTvTylbrAN0QHAQBUCkVYQV9GL5eaBkoPNBQdBABQKRRhRRiNxs/nRfJN9ABgAyjCiti6dWu2Uy2pUVvRQQAAlUURVsSM2eHZXSaKTgEAMAOKsNwSExP/vHxVeqqX6CAAADOgCMtt+udf5AWOZhE9ANgGirB8kpKSfv7pJ2OHoaKDAADMgyIsnznhEcb2r7CIHgBsBkVYDllZWdHLV+g6jRIdBABgNhRhOUQvX6Fp9rTk+ZDoIAAAszFDEZ47d653794PPvigj49P5femWAaD4fPwqOygsaKDAADMyQxFqNVqe/fu/emnn2ZmZlZ+b4q1devWXNc60sNtRAcBAJiTGdYANGrUqFGjRomJiZXflZLNmB2eFfSm6BQAADOz6GK4wsLCmzdvXrx40fTQycnJWt5NPXr06IUrSdLwENFBAABmZtEi/OOPP3788ceFCxeaHjo7O+/fv9/JycmSGSpm2szP8jqPlezsRQcBAJWTc3JysrOzy7i1s7OzVltK01m0CB9//PGxY8eGhoZa8qCVd/ny5bi4OOPMpaKDAAA0rq6ubm5uZtwjyydK90VElNxxmORSXXQQAID5mWEi1Ol0CQkJf/75p8Fg2LVrl4uLS4cOHSq/W4XIyspasXKVbvIh0UEAAFXCDEWYk5OzbNkySZJCQkKWLVvm5eVlS0W47MtozX+6SrV9RQcBAFQJMxRh7dq1N27cWPn9KJDBYJg9f0H2y6tFBwEAVBWuEZZky5Ytee71pIatRQcBAFQVirAk01lEDwC2jiK8ryNHjvx99br0RE/RQQAAVYgivK+PZs3J6zyORfQAYNsowuL9s4i+g5Wt/QcAlBdFWLw54RHGDsMkJ3N+eAEAQIEowmJkZWWtXLWab6IHADWgCIuxZOkyzePdpdr1RQcBAFQ5ivBuBoNhTsTCnE6jRQcBAFgCRXi3zZs362r6Sg39RQcBAFgCRXi3GXPm3w5kET0AqAVF+D8OHDhw+Xqa9EQP0UEAABZCEf6PmXPCc4PGSBpOCwCoBa/4/7p06VJ8fLyx3cuigwAALIci/NfnX8w3dBzOInoAUBWK8B+3b99evWZtQafXRQcBAFgURfiPRUuWalo8J9XyER0EAGBRFKEkSVJhYeHciAU5nd4QHQQAYGkUoSRJ0qZNm/QeD0sPtRIdBABgaRShJEnSjLkRLKIHAHWiCKX9+/f/90aG1OI50UEAACUyGox6nb29mb8vnSKUpn8+LydoLIvoAUDpjm199D/NPTw8zLtXtb/6//333/v375fbDhEdBABQCvf4iGnvTjD7btVehJ/NDTcEjJCcXEUHAQCU6NJRl5yUXr16mX3Hqi7CzMzMNWvXFgSMFB0EAFAK193zJk8cb/YLhJIkac2+RyuyaMlSTctgFtEDgNKlXZHO7R0xfGVV7Fu9E6Fer58XtSg3aKzoIACAUjjGRY149VU3tyr5LGj1ToQbN27UezaWfJ8QHQQAUKL8LPtf17617FgV7V69E+HMLyJvB44TnQIAUAq7hOhu3bs3aNCgivav0okwPj7+vzczpebPig4CACiR0eCSsGTK1m+q7ggqnQhnzg7P7fwmi+gBQOkStzR5yLd169ZVdwQ1NsGff/65/8ABuc0g0UEAAKVwj4+sikX0RamxCD+fF2EIDGMRPQAo3YVfXHJvBAcHV+lBVFeEGRkZ67/+Wt+JRfQAoHSueyPff+etqlhEX5TqbpZZuHiJXctgqUZd0UEAACVKuyyd3/fqsFVVfRx1TYR6vT58weKcQBbRA4DSOe2JDHttRBUtoi9KXRPhN998o/dqKvm2EB0EAFCivNt2h9aPjz5ugUOpayL8ZF5UVhCL6AFA6ewSvnz2ued8fX0tcCwVTYRxcXFXU29J/+kmOggAoETGQpf4xZO3b7bM0VQ0Ec6cHZ7bZQKL6AFA6RK3NH3kYX9/f8scTS2t8Oeffx48fFhuyyJ6AFA69/jID6t4EX1RainCz+aGFwa8Jjm4iA4CACjRhQNuuoyePXta7ICqKML09PSvN2zQ8030AKB4rnGRU94eb2dnuXpSxc0yUQsXa1qGSDW8RQcBAJQo9ZLmwv5hQ9dY8pi2PxHq9fqIRUtyA8eIDgIAKIVTXOTrI19zdbXoZ0Hb/kS4fv36wjrNWEQPAEqXm6k5uO7N5SctfFjbnwhnzJnPInoAUD77hOgewcH169e38HFtfCLcvXt3SnYBi+gBQOmMhc77lry/c6vlj2zjE+HMOeE5QeMkjUZ0EABAiY5sbta08ZNPPmn5I9tyEZ4/f/7QkaNym5dEBwEAlKL6vihLLqIvypaL8NM58/QBI1lEDwBKdz7BTX/r+eefF3Jwmy3C9PT0jRs3FQa8JjoIAKAUbnsjP3hngiUX0RdlszfLREQtlJ7qzSJ6AFC61EvSX7+88so6Uce3zYlQp9NFLFycFzhadBAAQCmc9swf/XqYhRfRF2WbE+H69euN9VtI9ZuLDgIAKFFuht2hb95cdVpgBPNMhCtXrnz++ef79eu3b98+s+ywkmbOjcgKZBE9ACid/b4ve4aE1K1bV2AGM0yE69evnzZtWnR0dHJycnBw8NGjRx955JHK77bCfv755xu5hdJjXQVmAACUzqB3Tlg65cdtYlOYoQgjIiI+/vjj7t27S5KUkJCwdOnSuXPnVn63FTZjdnhO4JssogcApTuy6fFmTVu2bCk2RWXfGpVl+cSJE+3atTM9bNeuXWJiYqVTVdy5c+eOJh6T2wwUmAEAUBbVExZME7SIvqjKToTp6ekFBQW1a9c2PfTw8EhJSbnfxmfPno2Li4uIiDA9dHNzi4mJcXJyqmSGoj765PPCwDDJwdmM+wQAmN+5eDf97Y4dO2ZnZ1fdQZydnbXaUpquskXo5uam0Wjy8vJMD3Nyctzd3e+3ccOGDZ9//vmePXuaHjo5OXl4eFQyQFE3b97cujVG/9EpM+4TAFAV3PZGTps0sXr16qKDVLoInZycvLy8Ll68+NBDD0mS9Pfffzdo0OB+Gzs7Oz/00EN+fn6VPOj9RC5YZOf/glS9ThXtHwBgHjcu2P19+JVXNojOIUlmWT4xcODAZcuWSZKUnZ29fv36gQPFXJ/T6XRRi5fmBo4VcnQAQNk57YkcPSrMxUURnwVthiKcMmXK+fPnH3vssSZNmvj7+/fu3bvy+6yAtWvXGn1bSvWaCTk6AKCscjPsjmwcN3qU6Bz/MMPyCS8vr8TExPPnz7u5ufn4+FR+hxXzyReRWc99LuroAIAyso9f2qt3b29vpXwWtHk+Yk2j0TRt2tQsu6qYn3766WauQXq0i8AMAIDSGfTO+79876dY0Tn+ZSMfuj1jdnhOlwksogcApTu8ocV/mj3xxBOic/zLForwzJkziceOy/4DRAcBAJTCPT5SCYvoi7KFIpw1d74+aBSL6AFA6f6IqyHlmT6SUzmsvghv3ry5ZcuWwk4jRQcBAJTCbW/ktElvaxR2Gcvqv48wcsFCyf8Fyc1TdBAAQIlS/rS7dHTIkE2ic9zNuidCnU4XtXhZXuAY0UEAAKVw2RMx9o3XFbKIvijrnghXr1ljbPCUVPdR0UEAACXKSZeObh6z7ozoHMWw7olw1ryorCC+iR4AlE4bv6Rv377KWURflBVPhDt37ryZZ5SaBokOAgAoUWGB476l7+7ZKTpH8ax4IpwxOzyny1ssogcApTv8zZMtn2jRooXoHMWz1iI8ffr0iVOnWUQPAMrnHh+ltEX0RVlrEc6aG17Q6XVJ6yg6CACgRGd3ezgYnnnmGdE57ssqi/DGjRsxMVsLA14THQQAUArXuMgP3p2gtEX0RVnlzTIRUQvl1gMkNw/RQQAAJbp+Xvvf44MHbRGdoyTWNxHqdLoFS5blB44WHQQAUArnPfPHjh7l7Kzoz4K2vonwq1Wr5If8JW+RX38IAChd1k3p6OYxX/8hOkcprGwilGX5k7kRLKIHAOXT7lv6Qv8BderUER2kFFY2Ef7www+ZRiepaaDoIACAEhXqHBO+nLxvl+gcpbOyiXDG7PDsoDdFpwAAlObXr/2eevKxxx4TnaN01lSEp0+fPnnmrNzqBdFBAAClcE9YoORF9EVZUxHO+PwLXeAoFtEDgNKd+dnDwfj000+LzlEmVlOEN27ciN22zdBxuOggAIBSuMVHTntvopIX0RdlNTfLzIuIktsMZBE9ACjd9fP2//3tpYFbRecoK+uYCPPy8hYtWZbf6Q3RQQAApXDZPW/82NEKX0RflHVMhF+tWiU1bid5NxEdBABQoqyb0vGtYzadE52jHKxgIpRl+dMvIrMCWUQPAEpnH79kQP8Bnp6eooOUgxVMhDt27LgluUhNAkQHAQCUqFDntD96UsJu0TnKxwomwumzw7OCxotOAQAohebgOv9Wfs2aNRMdpHyUXoSnTp06c/ac1Kqf6CAAgFK4JSy0lkX0RSm9CKd/NlcXNFqydxAdBABQotM/ejpJnTt3Fp2j3BRdhMnJydtjYw0dXxUdBABQCre9kR9OfttaFtEXpeginB+1UGo3WHKtLToIAKBESWe01868NHCg6BwVodwiNBqNi5fyTfQAYAVc9kaOH/OGo6NVfha0cpdP6PX6vJwsyauR6CAAgBJl3ZCOfzdm83nROSpIuRMhAMAqaOOXvDRwoIeHtX4WtHInQgCAFSjUOSZEv3MgTnSOimMiBABUnOaXNW3atH700UdFB6k4JkIAQEXJslt81LRVi0TnqBQmQgBARZ3e+YCrQ1BQkOgclUIRAgAqyG1v5EdT3hGdorIoQgBAhSSdcUj548UBA0TnqCyKEABQEdX2hE98c4yVLqIviptlAADll3VD/m3761vmic5hBkyEAIByc4hbNHjQIOtdRF8UEyEAoJz0edr90W8f3Cc6h3kwEQIAykfzy5r27ds3bdpUdBDzYCIEAJSHLLvtWzhtzRLROcyGiRAAUB6nvveu4dKpUyfROcyGIgQAlIPb3sgPJ70lOoU5UYQAgDJLOu1w43z//v1F5zAnihAAUFbVdoe/PX6sDSyiL4qbZQAAZZN5Tf5te1hMuOgcZsZECAAoE4f4xUOGDLGNRfRFMRECAMqgIFd7YMW74QdE5zA/80yEOp3uyJEjcXFxZtkbAEBpNL+s7tihQ+PGjUUHMT8zTIQJCQldu3b18PC4detWTk5O5XcIAFAWWXbbt2jqumWic1QJM0yELVu2TElJiY2NrfyuAABKdHJHvVquAQEBonNUCTNMhO7u7pXfCQBAsdz3Rnw4ZaLoFFXFojfL3Lx5c/v27UlJSaaH7u7uI0eOtLMrfirV6/WSbMFwAIBiXT3lkHqhV69eer1edJRys7e3v1/L3FGmIjx8+PDQoUPvfT42NrZRo0ZlD1RQUGA0GtPT000Ps7OzCwoKHBwcit3YYDCUfc8AgCrisvuLCWNH29nZWePLcqktKJWxCJs3b75169Z7n/f19S1XIB8fn86dO4eGhpZlY41GI2nKtXsAgLllXtOc2jkmdqGzs7PoKFWlTEXo4uLSpEmTqo4CAFAah70LX3755Zo1a4oOUoXMcI0wMzNz0qRJqampBQUFYWFhnp6en3zySeV3CwAQrCBXe2DlOxEHReeoWmYoQgcHBz8/P0mSunfvLnETKQDYCs2BrwI7dSrXvSDWyAxF6OrqOnLkyMrvBwCgILLRNWHxB18vF52jyvGh2wCA4pzYXr+2e4cOHUTnqHIUIQCgGG7xER9Nflt0CkugCAEA97h8zDnzSr9+/UTnsASKEABwt2pxEe+MH6vVquKr+lTxLwkAKIfMZOn0j2E7FovOYSFMhACA/+EYt2BoaGiNGjVEB7EQJkIAQBG6HPtfVk2M+lV0DsthIgQA/MvuwMqgoKCHH35YdBDLYSIEAPw/2egSv3Dq5jWic1gUEyEA4P8d39bAq3a7du1E57AoihAA8A+3vREfT3lHdApLowgBAJIkSdKlRJespD59+ojOYWkUIQBAkiTJNW7+pAnjVLKIvijV/QsDAIqRkST/vmvE90tF5xCAiRAAIDnGRb06dKh6FtEXxUQIAKqXn2X3y6q3Fh4RnUMMJkIAUDu7A1893bVrw4YNRQcRg4kQANTNaKi2b9EH364VnUMYJkIAULcT2x6u7922bVvROYShCAFA1dz2Rkx7d4LoFCJRhACgYpcSXbKSe/fuLTqHSBQhAKiX657wyRPH29vbiw4iEkUIAGqVdkX+fdfwV4eJziEYRQgAKuWwd8GI4cOrV68uOohgLJ8AAFXKz9IeXDNx6THROcRjIgQANbLbv6Jzly4NGjQQHUQ8JkIAUB+jodq+xe+sXy46hyIwEQKA+hzb2vhBn1atWonOoQgUIQCojnu82hfRF0URAoDKXDrqkpMSEhIiOodSUIQAoC6uu+dNeXuCyhfRF0URAoCapF2Rzu1lEX1RFCEAqIhjXNRrI4a7ubmJDqIgLJ8AANXIz7L/de1bXx4XnUNZmAgBQC3sEqK7P/usr6+v6CDKwkQIAOpgNFRLWDJ56zeicygOEyEAqEPilkcaNmjdurXoHIpDEQKAKrjHR37IIvriUIQAoAIXfnHTpQcHB4vOoUQUIQDYPte9kZMnvmlnx2t+MTgpAGDr0i5L5/cNGxoqOodCUYQAYOOc9kS+PvI1FtHfD8snAMCmZafZHVo/YcVvonMoFxMhANgyx13hA/r39/HxER1EuZgIAcB2ZadqD6yYcfKY6ByKxkQIADbL6cfZQwYP4jPVSsZECAA26tY1u4Orp0WfFJ1D6ShCALBN1WLeCwsbydXBUlGEAGCLftte+/qxTz5aLjqHFaAIAcDmZKdW+2bM1zEbXFxcREexAtwsAwC2xmXDm6+FDunYsaPoINaBiRAAbIrmRKxn6qlZM1aJDmI1KEIAsCHZqS4bxm7avoU3RcuOt0YBwHZU2zAu7NVX2rRpIzqINTHbRJifn6/VarVaRkwAEENzZJNn6plPP14jOoiVMcNEuHHjxmbNmtWuXbtGjRq9evVKS0ur/D4BAOWTnVptyzub169ydnYWHcXKmKEInZ2dV69enZOTk5KSkpeXN2nSpMrvEwBQLq7fjH3jtWH+/v6ig1gfM7yTGRISYvrBzc2tT58+a9eurfw+AQBlpzmy8YHMczM+XCc6iFUy5yU9WZZjYmI6d+58vw0MBsPNmzcvXrxoeujq6lqnTh0zBgAANcpMdtk8cfNPO5ycnERHsUoaWZZL3Sg1NXXhwoX3Ph8aGvrQQw/deThjxowNGzYcOnTI1dW12P34+flduXLlzrckOzk57d+//37/5XQ6nXe9+oZFWaXGAwD1ko2uUc+PDW4/5b13y/tHc3Jy7vdybTOcnZ0dHBxK3qasE6GdXSlXEyMiIlavXh0fH1/CaW3evPm4ceNCQ0PLckRHR0dJU8Z0AKBS2l3hTV31M6d/ZG9vX94/q9Fo7kwmalamIvT09Jw6dWoJG0RHR4eHh+/du7devXpmCgYAKM3lY867wr9NPFSBFsQdZrhGuHbt2jFjxkRFRV24cOHChQsuLi4dOnSo/G4BACXJu+26csjyJQuKXqJCBZihCFNSUgICAjZu3Gh66O3tTRECQNWSjdXWDB/Q45kB/fuLjmL1zFCEEydOnDhxYuX3AwAoI8fYj5o5ZCyJ2ig6iC3gE9EAwMpojm6udWLD90d/dXR0FJ3FFlCEAGBV/jpYbZxySjIAAAuzSURBVOObu/bt8fLyEh3FRvDtEwBgPdKuuHw58JvVKx5//HHRUWwHRQgAVkKX7bqkzydTJ/fs2VN0FJtCEQKANZCN1VaG9g5qPeHNsaKj2BquEQKAFXDc8l4zp1srlm4WHcQGUYQAoHR2v6z2+vOHn48c5DbRqkARAoCy/bnfLfaD3Qfia9WqJTqKbeIaIQAoWOqlaiuGbPlmbZMmTURHsVkUIQAoVd5t1yV9Pp8+9emnnxYdxZZRhACgSEaD61evvBLyzJhRr4uOYuMoQgBQIudNbz1VyxA5b47oILaPm2UAQHHs96/wurw39shBrZZX6SrHKQYAhfl9l+v2D3f/ur9GjRqio6gCRQgASnL9fLWvhn63ZUPjxo1FR1ELrhECgGLkpFdb0nveZzODgoJER1ERJkIAUAaD3jX6pZGDXggb+ZroKOrCRAgAiuC8Ybx/PZc5s2aKDqI6TIQAIJ79z+H1bxyJ/TXB3t5edBbVoQgBQLTTP7rHR+45ctDNzU10FDWiCAFAqOSzrmtG/LBjq6+vr+goKsU1QgAQJzvVdWm/pZHz2rZtKzqKejERAoAghQVu0QPHDR88ePAg0VFUjYkQAMRwWRfW6RGvmR9/KDqI2jERAoAA2h8+a3D77KYfEzQajegsakcRAoClaY5/V/2XZbuOHKxWrZroLKAIAcDCrpyo9vUbO3/cUb9+fdFRIElcIwQAi7p1rdrSF5YvWeDv7y86Cv7BRAgAlqLPc13Wf9Kbo14cMEB0FPyLIgSAKibL0qWjDidiHI9927Nr4NQp74kOhP9BEQJA1ZCN0pUTDqd2OB5ZX83O+ELv4EGTV3fo0EF0LNyNIgQAs5KN0l8HnY5/a38sxqNWjdCX+veeusHPz090LNyX4oswN0N0AgAoA1mWLiU6/RajOb6tvq9v6MB+A6L2NGnSRHQslE65RajVah/wqnN76qOig5SXLEksjy0jzlXZca7KTti5euiRpqEv9uu//NeGDRsKCYCKUW4R2tvbX/vvZdEpyi0rK8vd3V10CuvAuSo7zlXZca5QXqwjBACoGkUIAFA1ihAAoGoUoTmlpqYuX75cdArrUFBQMH/+fNEprMacOXNkWRadwjosWrQoOztbdArrsG7duqSkJNEpxKMIzenixYtr164VncI6ZGRkLFiwQHQKqzF79uyCggLRKaxDdHR0cnKy6BTWYePGjb///rvoFOJRhAAAVaMIAQCqRhECAFRNY8kr8G3btv3777/d3NwsdkQLKygoSE1NrVevnuggVsBoNF69erVBgwaig1iHS5cuPfjggxoNHy5TuqtXr3p7e2u1yv20EOW4fv16zZo1nZ2dRQepQoMGDZoxY0bJ21i0CDMzM1NTU+3sbHkM1el0Tk5OolNYB85V2XGuyo5zVXZqOFd169Z1cXEpeRuLFiEAAEpjy8MZAACloggBAKpGEQIAVI0iBACoGncYm8fFixcTExMzMjIGDhxYvXr1YreJiYlJSEioX7/+yJEjbXgNSalkWd6wYcPhw4cbNmw4YsSIe2/oOnDgwJkzZ+48HDFihG3faXyXQ4cObd682d3dfdiwYb6+vvducO3atZUrV6alpfXu3TsgIMDyCZXjypUrK1euzM7O7t+/f+vWre/6pwUFBV999dWdhy1btrx3G5XIyck5fvz4H3/80bhx46CgoGK3SU5OXrFiRWZmZp8+fTp06GDZgIKp6PWl6iQnJ7dq1WrJkiVhYWGpqanFbjN37tx33nmncePGCQkJXbt2VfPNulOnTp05c+YjjzyyY8eO3r1737vBhg0b1qxZc/H/WT6hQLt37+7evbu3t3d6erq/v//Nmzfv2uD27dtt2rS5fPmyr69vnz59tm3bJiSnEqSkpPj7+2dkZHh7e3fr1m3v3r13bZCbmxsWFvbXX3+ZfpHS09NFxFSEiRMnvv7667Nnz161alWxG2RkZLRu3To5OdnHxyckJOSHH36wcELBZFSawWCQZVmn00mS9Ndff927gU6nq1OnTnx8vCzLer3e19f3559/tnRKZbh9+3b16tV/++03WZbz8vJq1ap15MiRu7YZO3bsRx99JCKdeM8888wXX3xh+jk4OPjTTz+9a4MFCxZ06tTJ9HN0dHTbtm0tmk9JZsyYERISYvp57ty5zz777F0bZGRkSJJUWFho8WiKY3qNeu+994YOHVrsBuHh4U8//bTp58WLFwcEBFgunAIwEZpBqW/c/fHHH7dv3+7YsaMkSVqttnPnzvHx8RaJpjjHjx93cXFp0aKFJEnOzs4BAQHFnopjx47Nnj37m2++yc/Pt3hGYWRZTkhI6Natm+nhM888c+/JiY+PL7rBoUOHTH8DU6F9+/YVPRX79u0rdrMFCxZERUX99ttvFoymOKW+Rt31e3XgwIHCwsKqz6UUFKElXL9+3dPT887vYp06dVT7NTHXr19/4IEH7jws9lT4+PjUq1cvMzMzPDy8ZcuWmZmZls0oTEZGRn5+/p3z4+Xlde3atbu2uXbtWtENZFm+dxuVuOtU5Obm3vWrotFounXrduPGjdOnTwcEBERFRYmIaR3uOplGozElJUVsJEviZpkySUtLq1Onzr3Pf/311/379y/1j2u12qJ/vdLr9bb9sUYtWrS490vORo0aFRUVpdVqDQbDnSeLPRWTJk0y/WA0Gtu3b7948eLJkydXaWCFcHBwkCTpzq9KYWGho6PjvdsU3UCSpHu3UYmi/1uZfjCdwDtq1Kjx448/mn7u06dP3759w8LCVHu6Sqby3yuKsEw8PDwq80ZBvXr10tLS8vPzTR9um5SU1Lx5c/OlU5yTJ0/e7x/Vq1fv2rVrRqPRNB8nJSU99dRT99vYzs6uXbt26rlfxt3d3d3dPSkpycfHR5KkpKSkunXr3rWNj4/PnRn66tWrWq3Wy8vL0kGVoeipSEpKqlmzpqur6/027tChQ15e3vXr1/mc92LddTKdnJw8PDzERrIk3hqtQufOnTt79qwkSU2bNm3YsOF3330nSVJGRsauXbtCQkJEpxOjVatW1apV27VrlyRJ169fP3DgQM+ePSVJSklJOXjwoGmbO9cF8/Lyfv755//85z+i0lpecHDwpk2bJEkyGAxbtmwx/Z7o9fo9e/bk5eWZNti6dater5ckafPmzc8995xqv2YhODj422+/NRqNkiRt2rQpODjY9PyhQ4dMbxebzphJbGxszZo1TX/DgElBQcGePXtM15iDg4NjYmJMf93ftGlTjx49VLVmibtGzaNz586myebxxx/38/PLzc2VZXnkyJGhoaGmDWJiYjw9PYcOHdqsWbP73bilEqtXr/by8ho2bFjjxo3Hjx9venLdunUNGzY0/ezr69ujR4/Bgwf7+vp27tw5Ly9PXFhL+/333728vF588cWAgIDWrVubfpFMV2suXLggy3JBQUGnTp3atGkzePBgT0/PxMRE0ZGFycnJ8fPz69Sp04svvlinTp2zZ8+ann/00UeXL18uy3JkZGTz5s2HDBnSrVu36tWrb9y4UWhekdatW+fn5+ft7e3h4eHn5xcZGSnL8tWrVyVJunz5sizL+fn57du3b9eu3aBBgx544IETJ06IjmxRfPuEeZw4caLopa8nn3zSzs7u0qVLsiw3bNjQ9OSlS5cOHjzo6+trun1Uzc6fP5+YmPjwww+3adPG9Ex6evrVq1dNd5NevXo1MTExPz+/UaNGrVq1EppUgPT09F27dlWvXr1Lly6m6zSFhYXHjh174oknTNdTCwsL4+LiMjIygoKCVPu+qIlOp9uzZ092dnbXrl1r1aplevLUqVN169b19PQsKCg4evTo5cuXa9as6e/v7+npKTatQCkpKabaM/H29vbx8dHr9cePH2/ZsqXp10yv18fFxWVmZnbp0kVt54oiBACompreBQYA4B4UIQBA1ShCAICqUYQAAFWjCAEAqkYRAgBUjSIEAKgaRQgAUDWKEACgahQhAEDVKEIAgKr9H2sdxnnM46RKAAAAAElFTkSuQmCC",
"text/html": [
"\n",
"\n"
],
"image/svg+xml": [
"\n",
"\n"
]
},
"metadata": {},
"execution_count": 10
}
],
"cell_type": "code",
"source": [
"using Plots\n",
"plot((1 - α)^(-1) * Fs(10, 0))"
],
"metadata": {},
"execution_count": 10
},
{
"cell_type": "markdown",
"source": [
"---\n",
"\n",
"*This notebook was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*"
],
"metadata": {}
}
],
"nbformat_minor": 3,
"metadata": {
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.10.2"
},
"kernelspec": {
"name": "julia-1.10",
"display_name": "Julia 1.10.2",
"language": "julia"
}
},
"nbformat": 4
}