Lagrangian Approximations for the Stochastic Reachability of a Target Tube

This example will demonstrate the use of SReachTools for Lagrangian-based verification of stochastic continuous-state discrete-time linear time-invariant (LTI) systems. This example script is part of the SReachTools toolbox, which is licensed under GPL v3 or (at your option) any later version. A copy of this license is given in https://sreachtools.github.io/license/.

Contents

Lagrangian Methods

Lagrangian methods perform computations with sets using operations like unions, intersection, Minkowski addition/differences, etc. This computation using set operations can be used to approximate (either over or under) the stochastic reachability of a target tube problem. We will demonstrate that this approach, while being be approximative, can outperform the current state-of-the-art dynamic programming solution in terms of computation time.

Advantages:

Disadvantages:

The theory for this approach can be found in

Further, we explore multiple implementations of the Lagrangian-based verification, where the vertex-facet enumeration is mitigated either via recursion-free support method in overapproximation or vertex-complexity preserving support vector method in underapproximation.

All computations were performed using MATLAB on an Ubuntu OS running on a laptop with Intel i7 CPU with 2.1GHz clock rate and 8 GB RAM. For sake of clarity, all commands were asked to be verbose (via `SReachSetOptions`). In practice, this can be turned off.

% Prescript running: Initializing srtinit, if it already hasn't been initialized
close all;clearvars;srtinit;

Problem Definition

In this example we will look at the viability problem for a double integrator. The system dynamics are:

$$  x_{k+1} = \left[ \begin{array}{cc}    1 & T \\   0 & 1  \end{array}\right]
x_{k} + \left[\begin{array}{c}    \frac{T^{2}}{2} \\    T  \end{array}\right]
u_{k} + w_{k}$$

where $x_{k+1} \in \mathbf{R}^{2}$ is the state, $u_{k} \in \mathbf{R}$ is the input, and $w_{k} \in \mathbf{R}^{2}$ is the disturbance. The following code defines this system with $w_{k}$ as an i.i.d. Gaussian disturbance with mean $[0, 0]^{\top}$ and variance $\mathrm{diag}(0.01, 0.01)$.

example parameters

T = 0.25;

% define the system
sys = getChainOfIntegLtiSystem(2, ...
    T, ...
    Polyhedron('lb', -0.1, 'ub', 0.1), ...
     RandomVector('Gaussian', zeros(2,1), 0.001*eye(2)));

Viability problem as a stochastic reachability of a target tube problem

We examine the viability problem in which we are interested in staying in a set of safe states. In this example the safe set is $\{x \in \mathbf{R}^{2}: |x_{i}| < 1, i = 1, 2\}$. The stochastic reachability of a target tube problem posed as a viability problem by constructing a target tube in which all sets in the tube are the safe set.

time_horizon = 5;

% safe set definition
safe_set = Polyhedron('lb', [-1, -1], 'ub', [1, 1]);
% target tube definition
target_tube = Tube('viability', safe_set, time_horizon);
% probability threshold desired
beta = 0.8;
% Plotting of target tube
figure(1)
clf
hold on
for time_indx = 0:time_horizon
    target_tube_at_time_indx = Polyhedron('H',...
        [target_tube(time_indx+1).A, ...
        zeros(size(target_tube(time_indx+1).A,1),1), ...
        target_tube(time_indx+1).b], 'He',[0 0 1 time_indx]);
    plot(target_tube_at_time_indx, 'alpha',0.25);
end
axis([-1 1 -1 1 0 time_horizon]);
box on;
grid on;
xlabel('x');
ylabel('y');
zlabel('time');
title('Target tube');
axis equal;

Lagrangian underapproximation for stochastic reachability of a target tube

SReachSet can compute Lagrangian underapproximation via multiple approaches. It can approximate the disturbance set, against which the robust computation is done, as a polytope or an ellipsoid (as specified by bound_set_method). Further, it can either perform an exact (and time-consuming) vertex-facet enumeration or an underapproximative (and complexity-preserving) computation using support vectors (as specified by compute_style). The vertex-facet enumeration may be performed by CDD or LRS, two popular enumeration techinques. We use CDD by default, but it can be switched to LRS using vf_enum_method.

n_dim = sys.state_dim + sys.input_dim; % Require unit vectors to sample X x U
lagrange_under_time = zeros(4,1);

Underapprox. type 1: Bound_set_method - Ellipsoid | Compute_style - VFmethod

This is the best method in this case, since it leverages the fact that the distubance is Gaussian and for two-dimensional polytopes, vertex-facet enumeration is not significantly hard.

timerVal=tic;
luOpts1 = SReachSetOptions('term', 'lag-under', 'bound_set_method', ...
    'ellipsoid', 'verbose',1,'compute_style','vfmethod');
luOpts_time(1) = toc(timerVal);
timerVal=tic;
luSet(1) = SReachSet('term', 'lag-under', sys, beta, target_tube, luOpts1);
lagrange_under_time(1) = toc(timerVal);
Computing Lagragian under approximation

Time_horizon: 5
Computation for time step: 4
Computation for time step: 3
Computation for time step: 2
Computation for time step: 1
Computation for time step: 0

Underapprox. type 2: Bound_set_method - Ellipsoid | Compute_style - Support

Compared to type 1, this method will result in more conservative solution. The benefit of this approach becomes clear when we use it in high-dimensional systems where vertex-facet enumeration is hard. By the virtue of being recursion-free, this approach will scale better than type 1.

timerVal=tic;
luOpts2 = SReachSetOptions('term', 'lag-under', 'bound_set_method', ...
    'ellipsoid', 'system', sys, 'n_vertices', 2^n_dim*10+2*n_dim,...
    'verbose',1,'compute_style','support');
luOpts_time(2) = toc(timerVal);
timerVal=tic;
luSet(2) = SReachSet('term', 'lag-under', sys, beta, target_tube, luOpts2);
lagrange_under_time(2) = toc(timerVal);
Spreading 86 unit-length vectors in 3-dim space
Analyzing 10 unit-length vectors in first quadrant
 1. Setting up the CVX problem...
  1 |  2 |  3 |  4 |  5 |  6 |  7 |  8 |  9 | 10 |
Solving the CVX problem...done
Status: Solved
Sum of slack: 1.568e-09 (< 1.000e-08)
Change in opt cost: 3.010e-01 (< 1.000e-05)

 2. Setting up the CVX problem...
  1 |  2 |  3 |  4 |  5 |  6 |  7 |  8 |  9 | 10 |
Solving the CVX problem...done
Status: Solved
Sum of slack: 3.182e-10 (< 1.000e-08)
Change in opt cost: 4.522e-02 (< 1.000e-05)

 3. Setting up the CVX problem...
  1 |  2 |  3 |  4 |  5 |  6 |  7 |  8 |  9 | 10 |
Solving the CVX problem...done
Status: Solved
Sum of slack: 1.659e-11 (< 1.000e-08)
Change in opt cost: 1.056e-02 (< 1.000e-05)

 4. Setting up the CVX problem...
  1 |  2 |  3 |  4 |  5 |  6 |  7 |  8 |  9 | 10 |
Solving the CVX problem...done
Status: Solved
Sum of slack: 3.766e-11 (< 1.000e-08)
Change in opt cost: 5.915e-03 (< 1.000e-05)

 5. Setting up the CVX problem...
  1 |  2 |  3 |  4 |  5 |  6 |  7 |  8 |  9 | 10 |
Solving the CVX problem...done
Status: Solved
Sum of slack: 3.092e-10 (< 1.000e-08)
Change in opt cost: 9.338e-04 (< 1.000e-05)

 6. Setting up the CVX problem...
  1 |  2 |  3 |  4 |  5 |  6 |  7 |  8 |  9 | 10 |
Solving the CVX problem...done
Status: Solved
Sum of slack: 1.669e-11 (< 1.000e-08)
Change in opt cost: 3.288e-04 (< 1.000e-05)

 7. Setting up the CVX problem...
  1 |  2 |  3 |  4 |  5 |  6 |  7 |  8 |  9 | 10 |
Solving the CVX problem...done
Status: Solved
Sum of slack: 5.145e-11 (< 1.000e-08)
Change in opt cost: 9.588e-08 (< 1.000e-05)

Completed spreading the vectors!
Computing Lagragian under approximation

Time_horizon: 5
Computation for time step: 4
Computation for time step: 3
Computation for time step: 2
Computation for time step: 1
Computation for time step: 0

Underapprox. type 3: Bound_set_method - Polytope | Compute_style - VFmethod

Compared to type 1, this method will result in more conservative solution since the bounded disturbance set is larger in volume. This approach is best used when the disturbance is non-Gaussian. This approach is best used when the disturbance is non-Gaussian.

luOpts3 = SReachSetOptions('term', 'lag-under', 'bound_set_method', ...
    'polytope', 'verbose', 1, 'compute_style','vfmethod','template_polytope',...
    Polyhedron('lb',-ones(sys.dist.dim,1),'ub',ones(sys.dist.dim,1)));
luOpts_time(3) = toc(timerVal);
timerVal=tic;
luSet(3) = SReachSet('term', 'lag-under', sys, beta, target_tube, luOpts3);
lagrange_under_time(3) = toc(timerVal);
Computing Lagragian under approximation

Computing the bounded disturbance set
Time_horizon: 5
Computation for time step: 4
Computation for time step: 3
Computation for time step: 2
Computation for time step: 1
Computation for time step: 0

Underapprox. type 3: Bound_set_method - Polytope | Compute_style - Support

Compared to type 1, this method will result in more conservative solution since the bounded disturbance set is larger in volume. Compared to type 3, this set will be more conservative since an underapproximative vertex-facet enumeration is performed.

luOpts4 = SReachSetOptions('term', 'lag-under', 'bound_set_method', ...
    'polytope', 'system', sys, 'n_vertices', 2^n_dim*10+2*n_dim,...
    'verbose', 1, 'compute_style','support', 'template_polytope',...
    Polyhedron('lb',-ones(sys.dist.dim,1),'ub',ones(sys.dist.dim,1)));
luOpts_time(4) = toc(timerVal);
timerVal=tic;
luSet(4) = SReachSet('term', 'lag-under', sys, beta, target_tube, luOpts4);
lagrange_under_time(4) = toc(timerVal);
Spreading 86 unit-length vectors in 3-dim space
Analyzing 10 unit-length vectors in first quadrant
 1. Setting up the CVX problem...
  1 |  2 |  3 |  4 |  5 |  6 |  7 |  8 |  9 | 10 |
Solving the CVX problem...done
Status: Solved
Sum of slack: 2.621e-10 (< 1.000e-08)
Change in opt cost: 2.621e-01 (< 1.000e-05)

 2. Setting up the CVX problem...
  1 |  2 |  3 |  4 |  5 |  6 |  7 |  8 |  9 | 10 |
Solving the CVX problem...done
Status: Solved
Sum of slack: 7.852e-10 (< 1.000e-08)
Change in opt cost: 7.511e-02 (< 1.000e-05)

 3. Setting up the CVX problem...
  1 |  2 |  3 |  4 |  5 |  6 |  7 |  8 |  9 | 10 |
Solving the CVX problem...done
Status: Solved
Sum of slack: 3.660e-10 (< 1.000e-08)
Change in opt cost: 1.453e-02 (< 1.000e-05)

 4. Setting up the CVX problem...
  1 |  2 |  3 |  4 |  5 |  6 |  7 |  8 |  9 | 10 |
Solving the CVX problem...done
Status: Solved
Sum of slack: 1.079e-10 (< 1.000e-08)
Change in opt cost: 8.359e-03 (< 1.000e-05)

 5. Setting up the CVX problem...
  1 |  2 |  3 |  4 |  5 |  6 |  7 |  8 |  9 | 10 |
Solving the CVX problem...done
Status: Solved
Sum of slack: 3.053e-11 (< 1.000e-08)
Change in opt cost: 4.323e-04 (< 1.000e-05)

 6. Setting up the CVX problem...
  1 |  2 |  3 |  4 |  5 |  6 |  7 |  8 |  9 | 10 |
Solving the CVX problem...done
Status: Solved
Sum of slack: 4.489e-11 (< 1.000e-08)
Change in opt cost: 7.156e-07 (< 1.000e-05)

Completed spreading the vectors!
Computing Lagragian under approximation

Computing the bounded disturbance set
Time_horizon: 5
Computation for time step: 4
Computation for time step: 3
Computation for time step: 2
Computation for time step: 1
Computation for time step: 0

Lagrangian overapproximation for stochastic reachability of a target tube

SReachSet can compute Lagrangian overapproximation via multiple approaches. It can approximate the disturbance set which augments the input space as a polytope or an ellipsoid (as specified by bound_set_method). Further, it can perform a recursive computation using vertex-facet enumeration or a recursion-free computation using support functions (as specified by compute_style). The vertex-facet enumeration may be performed by CDD or LRS, two popular enumeration techinques. We use CDD by default, but it can be switched to LRS using vf_enum_method.

n_dim_over = sys.state_dim;     % Require unit vectors to sample X
lagrange_over_time = zeros(4,1);

Overapprox. type 1: Bound_set_method - Ellipsoid | Compute_style - VFmethod

This is the best method in this case, since it leverages the fact that the distubance is Gaussian and for two-dimensional polytopes, vertex-facet enumeration is not significantly hard.

timerVal=tic;
loOpts1 = SReachSetOptions('term', 'lag-over', 'bound_set_method', ...
    'ellipsoid', 'verbose', 1, 'compute_style','vfmethod');
loOpts_time(1) = toc(timerVal);
timerVal=tic;
loSet(1) = SReachSet('term', 'lag-over', sys, beta, target_tube, loOpts1);
lagrange_over_time(1) = toc(timerVal);
Computing Lagragian over approximation

Time_horizon: 5
Computation for time step: 4
Computation for time step: 3
Computation for time step: 2
Computation for time step: 1
Computation for time step: 0

Overapprox. type 2: Bound_set_method - Ellipsoid | Compute_style - Support

Compared to type 1, this method will provide an overapproximative solution but can potentially be faster since it is recursion-free.

timerVal=tic;
loOpts2 = SReachSetOptions('term', 'lag-over', 'bound_set_method', ...
    'ellipsoid', 'verbose', 1, 'compute_style','support', 'system', sys,...
    'n_vertices', 2^n_dim_over * 7+2*n_dim_over);
loOpts_time(2) = toc(timerVal);
timerVal=tic;
loSet(2) = SReachSet('term', 'lag-over', sys, beta, target_tube, loOpts2);
lagrange_over_time(2) = toc(timerVal);
Spreading 32 unit-length vectors in 2-dim space
Analyzing 7 unit-length vectors in first quadrant
 1. Setting up the CVX problem...
  1 |  2 |  3 |  4 |  5 |  6 |  7 |
Solving the CVX problem...done
Status: Solved
Sum of slack: 2.188e-02 (< 1.000e-08)
Change in opt cost: 1.561e-01 (< 1.000e-05)

 2. Setting up the CVX problem...
  1 |  2 |  3 |  4 |  5 |  6 |  7 |
Solving the CVX problem...done
Status: Solved
Sum of slack: 5.575e-11 (< 1.000e-08)
Change in opt cost: 3.901e-02 (< 1.000e-05)

 3. Setting up the CVX problem...
  1 |  2 |  3 |  4 |  5 |  6 |  7 |
Solving the CVX problem...done
Status: Solved
Sum of slack: 3.749e-10 (< 1.000e-08)
Change in opt cost: 3.079e-03 (< 1.000e-05)

 4. Setting up the CVX problem...
  1 |  2 |  3 |  4 |  5 |  6 |  7 |
Solving the CVX problem...done
Status: Solved
Sum of slack: 1.058e-11 (< 1.000e-08)
Change in opt cost: 4.595e-05 (< 1.000e-05)

 5. Setting up the CVX problem...
  1 |  2 |  3 |  4 |  5 |  6 |  7 |
Solving the CVX problem...done
Status: Solved
Sum of slack: 1.848e-11 (< 1.000e-08)
Change in opt cost: 8.840e-09 (< 1.000e-05)

Completed spreading the vectors!
Computing Lagragian over approximation

Evaluating support function:    32/   32

Overapprox. type 3: Bound_set_method - Polytope | Compute_style - VFmethod

Compared to type 1, this method will result in more conservative solution since the bounded disturbance set is larger in volume.

timerVal=tic;
loOpts3 = SReachSetOptions('term', 'lag-over', 'bound_set_method', ...
    'polytope', 'verbose', 1, 'template_polytope',...
    Polyhedron('lb',-ones(sys.dist.dim,1),'ub',ones(sys.dist.dim,1)),...
    'compute_style','vfmethod');
loOpts_time(3) = toc(timerVal);
timerVal=tic;
loSet(3) = SReachSet('term', 'lag-over', sys, beta, target_tube, loOpts3);
lagrange_over_time(3) = toc(timerVal);
Computing Lagragian over approximation

Computing the bounded disturbance set
Time_horizon: 5
Computation for time step: 4
Computation for time step: 3
Computation for time step: 2
Computation for time step: 1
Computation for time step: 0

Overapprox. type 4: Bound_set_method - Polytope | Compute_style - Support

Compared to type 1, this method will result in more conservative solution since the bounded disturbance set is larger in volume. Compared to type 3, this method will provide an overapproximative solution but can potentially be faster since it is recursion-free.

loOpts4 = SReachSetOptions('term', 'lag-over', 'bound_set_method', ...
    'polytope', 'verbose', 1, 'template_polytope',...
    Polyhedron('lb',-ones(sys.dist.dim,1),'ub',ones(sys.dist.dim,1)),...
    'compute_style','support', 'system', sys,...
    'n_vertices', 2^n_dim_over* 7 +2*n_dim_over);
loOpts_time(4) = toc(timerVal);
timerVal=tic;
loSet(4) = SReachSet('term', 'lag-over', sys, beta, target_tube, loOpts4);
lagrange_over_time(4) = toc(timerVal);
Spreading 32 unit-length vectors in 2-dim space
Analyzing 7 unit-length vectors in first quadrant
 1. Setting up the CVX problem...
  1 |  2 |  3 |  4 |  5 |  6 |  7 |
Solving the CVX problem...done
Status: Solved
Sum of slack: 2.242e-02 (< 1.000e-08)
Change in opt cost: 1.447e-01 (< 1.000e-05)

 2. Setting up the CVX problem...
  1 |  2 |  3 |  4 |  5 |  6 |  7 |
Solving the CVX problem...done
Status: Solved
Sum of slack: 4.952e-09 (< 1.000e-08)
Change in opt cost: 4.312e-02 (< 1.000e-05)

 3. Setting up the CVX problem...
  1 |  2 |  3 |  4 |  5 |  6 |  7 |
Solving the CVX problem...done
Status: Solved
Sum of slack: 1.044e-10 (< 1.000e-08)
Change in opt cost: 9.922e-03 (< 1.000e-05)

 4. Setting up the CVX problem...
  1 |  2 |  3 |  4 |  5 |  6 |  7 |
Solving the CVX problem...done
Status: Solved
Sum of slack: 6.872e-10 (< 1.000e-08)
Change in opt cost: 5.405e-04 (< 1.000e-05)

 5. Setting up the CVX problem...
  1 |  2 |  3 |  4 |  5 |  6 |  7 |
Solving the CVX problem...done
Status: Solved
Sum of slack: 4.490e-11 (< 1.000e-08)
Change in opt cost: 1.555e-06 (< 1.000e-05)

Completed spreading the vectors!
Computing Lagragian over approximation

Computing the bounded disturbance set
Evaluating support function:    32/   32

Dynamic programming solution

We compare the results with dynamic programming to see how the approximations appear and how they compare in simulation times.

dyn_prog_xinc = 0.025;
dyn_prog_uinc = 0.1;
tic;
[prob_x, cell_of_xvec] = SReachDynProg('term', sys, dyn_prog_xinc, ...
    dyn_prog_uinc, target_tube);
dynprog_time = toc();

Compute the beta-stochastic level set

dyn_soln_lvl_set=getDynProgLevelSets2D(cell_of_xvec, prob_x, beta, target_tube);

Simulation times: Lagrangian approximation beats dynamic programming

The simulation times for Lagrangian computation is much faster (in most cases) than dynamic programming. Further, dynamic programming solution provides no approximation guarantees while Lagrangian approach provides grid-free approximation guarantee.

fprintf('Simulation times [seconds]:\n');
fprintf(' Lagrangian:\n');
% fprintf('   Overapproximation  : %.3f (online: %1.3f | offline: %1.3f)\n',...
%     lagrange_over_time + loOpts_time, lagrange_over_time, loOpts_time);
fprintf('   Overapproximation   online  |  offline | Total\n');
fprintf('(Ellipsoid,VFmethod)  %1.2e | %1.2e | %1.2f\n',...
    lagrange_over_time(1),loOpts_time(1),lagrange_over_time(1)+loOpts_time(1));
fprintf('(Ellipsoid, support)  %1.2e | %1.2e | %1.2f\n',...
    lagrange_over_time(2),loOpts_time(2),lagrange_over_time(2)+loOpts_time(2));
fprintf('(Polytope ,VFmethod)  %1.2e | %1.2e | %1.2f\n',...
    lagrange_over_time(3),loOpts_time(3),lagrange_over_time(3)+loOpts_time(3));
fprintf('(Polytope , support)  %1.2e | %1.2e | %1.2f\n',...
    lagrange_over_time(4),loOpts_time(4),lagrange_over_time(4)+loOpts_time(4));
fprintf('   Underapproximation  online  |  offline | Total\n');
fprintf('(Ellipsoid,VFmethod)  %1.2e | %1.2e | %1.2f\n',...
    lagrange_under_time(1),luOpts_time(1),lagrange_under_time(1)+luOpts_time(1));
fprintf('(Ellipsoid, support)  %1.2e | %1.2e | %1.2f\n',...
    lagrange_under_time(2),luOpts_time(2),lagrange_under_time(2)+luOpts_time(2));
fprintf('(Polytope ,VFmethod)  %1.2e | %1.2e | %1.2f\n',...
    lagrange_under_time(3),luOpts_time(3),lagrange_under_time(3)+luOpts_time(3));
fprintf('(Polytope , support)  %1.2e | %1.2e | %1.2f\n',...
    lagrange_under_time(4),luOpts_time(4),lagrange_under_time(4)+luOpts_time(4));
fprintf('   Dynamic programming: %.3f\n', dynprog_time);
Simulation times [seconds]:
 Lagrangian:
   Overapproximation   online  |  offline | Total
(Ellipsoid,VFmethod)  5.90e-01 | 2.41e-03 | 0.59
(Ellipsoid, support)  9.95e+00 | 8.97e+00 | 18.91
(Polytope ,VFmethod)  1.05e+00 | 8.62e-03 | 1.06
(Polytope , support)  7.15e+00 | 9.46e+00 | 16.62
   Underapproximation  online  |  offline | Total
(Ellipsoid,VFmethod)  4.29e-01 | 3.63e-03 | 0.43
(Ellipsoid, support)  2.84e+00 | 2.52e+01 | 28.00
(Polytope ,VFmethod)  9.86e-01 | 2.92e+00 | 3.91
(Polytope , support)  3.39e+00 | 2.15e+01 | 24.85
   Dynamic programming: 94.858

Plotting all the sets together

As expected, the over-approximation and the under-approximation obtained via Lagrangian approach bounds the dynamic programming solution from "inside" and "outside".

figure(2);
clf
plot(safe_set, 'color', 'k');
hold on;
plot(loSet(2), 'color', 'y','alpha',1);
plot(loSet(1), 'color', 'r','alpha',0.5);
plot(dyn_soln_lvl_set,'color', 'b')
plot(luSet(1), 'color', 'm','alpha',1);
plot(luSet(2), 'color', 'g','alpha',0.75);
hold off;
xlabel('$x_1$', 'Interpreter', 'latex')
ylabel('$x_2$', 'Interpreter', 'latex')
leg = legend('Safe set',...
    'Overapproximation (Ell, support)',...
    'Overapproximation (Ell, VF)',...
     'Dyn. prog. soln.',...
    'Underapproximation (Ell, VF)',...
    'Underapproximation (Ell, support)');
set(leg,'Location','EastOutside');
box on;
axis equal;
axis tight;

figure(3);
clf
plot(safe_set, 'color', 'k');
hold on;
plot(loSet(4), 'color', 'y','alpha',1);
plot(loSet(3), 'color', 'r','alpha',0.5);
plot(dyn_soln_lvl_set,'color', 'b')
plot(luSet(3), 'color', 'm','alpha',1);
plot(luSet(4), 'color', 'g','alpha',0.75);
hold off;
xlabel('$x_1$', 'Interpreter', 'latex')
ylabel('$x_2$', 'Interpreter', 'latex')
leg = legend('Safe set',...
    'Overapproximation (Poly, support)',...
    'Overapproximation (Poly, VF)',...
     'Dyn. prog. soln.',...
    'Underapproximation (Poly, VF)',...
    'Underapproximation (Poly, support)');
set(leg,'Location','EastOutside');
box on;
axis equal;
axis tight;