issue #54: Added code, added picture of the reachability tube and set

doc/chap_examples.rst

@@ -632,15 +632,13 @@ Defining :math:`x_1 = x, x_2=\dot{x}_1, x_3 = \theta` and :math:`x_4=\dot{x}_3`,
    0 \\
    -\frac{m_2Lg\pi}{4J_c+m_2L^2} \end{array}\right].   
 
-Consider some moment of time :math:`t_1` and final position :math:`x_1(t_1) = x_1, x_2(t_1) = 0, x_3(t_1) = \frac{\pi}{2}, x_4(t_1) = 0`. It is required to calculate the backward reachability sets (tube) for the linearized system :eq:`invpendls` emanating from the given final position and project it onto :math:`(x_1, x_3)` subspace. It is also required to identify whether it’s possible to reach the final position from a given initial position :math:`x_1(t_1) = x_1, x_2(t_1) = v_1, x_3(t_1) = \theta_1, x_4(t_1) = \omega_1` using some admissible control function.
-
-.. raw:: html
+Consider some moment of time :math:`t_1` and final position :math:`x_1(t_1) = x^1_1, x_2(t_1) = 0, x_3(t_1) = \frac{\pi}{2}, x_4(t_1) = 0`. It is required to calculate the backward reachability sets (tube) for the linearized system :eq:`invpendls` emanating from the given final position and project it onto :math:`(x_1, x_3)` subspace. It is also required to identify whether it’s possible to reach the final position from a given initial position :math:`x_1(t_0) = x^0_1, x_2(t_0) = v_1, x_3(t_0) = \theta_1, x_4(t_0) = \omega_1` using some admissible control function.
 
-	<h2>References</h2>
+.. literalinclude:: ../products/+elltool/+doc/+snip/s_chapter06_section05_snippet01.m
+   :language: matlab
+   :linenos:
 
-.. [SUN2003] L.Muñoz, X.Sun, R.Horowitz, and L.Alvarez. 2003. Traffic Density
-   Estimation with the Cell Transmission Model. In *Proceedings of the
-   American Control Conference*, 3750–3755. Denver, Colorado, USA.
+.. _brset:
 
 Pendulum with a flywheel
 -------------------------

products/+elltool/+doc/+snip/s_chapter06_section05_snippet01.m

@@ -0,0 +1,48 @@
+alpha = 10;
+k1 = 50;
+k2 = 500;
+m1 = 150;
+m2 = 100;
+L = 5;
+g = 9.8;
+Jc = 100;
+
+x1 = 1;
+v1 = 0;
+theta1 = pi/2;
+omega1 = 0;
+t1 = 1;
+
+d = Jc + m2*L^2/4;
+% define matrices aMat, bMat, and control bounds uBoundsEllObj:
+aMat = [0, 1, 0, 0; 0, -k1/m1, 0, 0; 0, 0, 0, 1; 0, 0, m2*L*g/2/d, -k2*L/d];
+bMat = [0; 1/m1; 0; 0];
+uBoundsEllObj = alpha * ell_unitball(1);
+%define disturbance:
+gMat = [0; 0; 0; -m2*L*g*pi/4/d];
+vEllObj = ellipsoid(1, 0); %known disturbance
+
+%linear system
+lsys = elltool.linsys.LinSysContinuous(aMat, bMat, uBoundsEllObj,...
+    gMat, vEllObj);
+timeVec = [t1, 0];
+% initial directions (some random vectors in R^4):
+dirsMat = [1 0 1 0; 1 -1 0 0; 0 -1 0 1; 1 1 -1 1; -1 1 1 0; -2 0 1 1].';
+% x1EllObj = ellipsoid([x1; v1; theta1; omega1], eps*eye(4)); %known final point
+x1EllObj = ellipsoid([x1; v1; theta1; omega1], zeros(4)); %known final point
+
+
+%backward reach set
+brsObj = elltool.reach.ReachContinuous(lsys, x1EllObj, dirsMat, timeVec,...
+    'isRegEnabled', true, 'isJustCheck', false, 'regTol', 1e-4,...
+    'absTol', 1e-5, 'relTol', 1e-4);
+
+basisMat = [1 0 0 0; 0 0 1 0].'; % orthogonal basis of (x1, x3) subspace
+psObj = brsObj.projection(basisMat); % reach set projection
+
+% plot projection of reach set internal approximation:
+psObj.plotByIa('b');
+
+%plot backward reach set approximation at time t=0
+psCutObj = psObj.cut(0);
+psCutObj.plotByIa('r');