Real Time Reachability Analysis for Marine Vessels

(1)

IN

DEGREE PROJECT INFORMATION AND COMMUNICATION TECHNOLOGY,

SECOND CYCLE, 30 CREDITS STOCKHOLM SWEDEN 2018,

Real Time Reachability

Analysis for Marine Vessels

SUDAKSHIN GANESAN

KTH ROYAL INSTITUTE OF TECHNOLOGY

SCHOOL OF ELECTRICAL ENGINEERING AND COMPUTER SCIENCE

(2)

(3)

i

Abstract

Safety verification of continuous dynamical systems require the computation of the reachable set. The reachable set comprises those states the system can reach at a specific point in time. The present work aims to compute this reachable set for the marine vessel, in the presence of uncertainties in the dynamic modeling of the system and in the presence of external disturbances in the form of wind, waves and currents. The reachable set can then be used to check if the vessel collides with an obstacle. The dynamic model used is that of a nonlinear maneuvering model for the marine vessel. The dynamics on the azipod actuators are also considered.

Several methods are considered to solve the reachability problem for the marine vessel. The first method considered is that of the Hamilton Jacobi Reachability analysis, where a dynamic game between the control input and the disturbance input is played. This results in a dynamic programming problem known as the Hamilton Jacobi Bellman Isaacs (HJBI) equation. It is solved using the Level-Set method, but it suffers from the curse of dimensionality. The other method considered is the use of set-theoretic approach, where an over-approximation of the reachable set is computed, in the context of safety verification. But on the downside, large sets of admissible control yields highly over-approximated reachable sets, which cannot be used

In order to overcome the disadvantages posed by the first two methods, emphasizing on the real-time computation, a third method is developed, where a supervised classification algorithm is used to compute the reachable set boundary. The dataset required for the classification algorithm is computed by solving a 2 Point Boundary Value Optimal Control Problem for the marine vessel. The features for classification algorithm can be extended, so as to include the uncertainties and disturbances in the system. The computation time is greatly reduced and the accuracy of the method is comparable to the exact reachable set computation.

Keywords: Reachable set, Safety Verification, Marine dynamics, Machine Learning, Optimal control.

(4)

ii

Sammanfattning

Säkerhetsverifiering av kontinuerliga dynamiska system kräver beräkning av mängden av tillstånd som kan nås vid en specifik tidpunkt, givet dess ini- tialtillstånd. Detta arbete fokuserar påatt bestämma denna mängd av nåbara tillstånd för ett marint fartyg under modellosäkerheter och externa störningar i form av vind, vågor och strömmar. Den nåbara mängden av tillstånd an- vänds sedan för att kontrollera om fartyget riskerar att kollidera med hinder.

Den dynamiska modell som används i våra studier är en icke-linjär modell där även dynamiken hos azipod-ställdonen betraktas.

Arbetet studerar flera metoder för att lösa problemet: en klassisk Hamilton- Jacobi nåbarhetsanalys, en mängd-teoretisk teknik, samt en ny metod baserad påmaskininlärning. Numeriska simuleringsstudier bekräftar att den föreslagna maskininlärningsmetoden är snabbare än de tvåalternativen.

Keywords: Säkerhetsverifiering, nåbarhetsanalys, fartyg dynamisk

(5)

iii

Acknowledgements

This thesis project was carried out at ABB Corporate Research Center, at Västerås, Sweden, at the Control, Optimization and Analytics group. It was supervised by the Division of Automatic Control at the Royal Institute of Technology (KTH) in Stockholm, Sweden.

I would like to thank my supervisors Hamid Reza Feyzmahdavian and Winston Garcia Gabin at ABB for giving me the freedom to explore the problem statement and the solution methods, which gave me a flavor for research in an industry and also for their continuous support throughout the project. I would like to thank my Supervisor and Examiner Mikael Johansson at KTH for giving me the opportunity to carry out this interesting project and also importantly for steering the research in an interesting direction.

At the end I would like to express my biggest gratitude to my loved ones, my family and friends, who have supported me throughout the entire process. I will be grateful forever for your love.

Sudakshin Ganesan,

Stockholm, 15th August 2018

(6)

List of Figures

1.1 Verification of the vessel’s safety using the Forward reachable set. . . 3 2.1 Exponential increase in complexity as the number of dimensions in the

state space increases [1]. . . 6 3.1 Motion in 6 degrees of freedom . . . 10 3.2 Illustration of the level set method, which encodes the target set and

the reachable set using value functions. . . 12 3.3 Forward Reachable set for the velocity states in the surge and the sway

direction. . . 15 3.4 Illustration of the overaproximation of reachable set used for safety ver-

ification [2]. . . 16 3.5 Illustration of the steps involved in approximate reachable set calculation

[2]. . . 16 3.6 Reachable set for the marine vessel with 3 DOF maneuvering model, for

3 min into the future, admissible control being 1/10^ththe original value. 17 4.1 Simplified 2D illustration of a cost limited reachable set for a given

threshold cost of J0 . . . 19 4.2 Illustration of the multiple shooting algorithm [3] . . . 21 4.3 Illustration of separating hyperplane & support vectors . . . 24 5.1 The Reachable set for the marine vessel with initial surge velocity, u

= 1 m/s. The boundary between the reachable and non reachable is classified using the SVM algorithm. . . 30 5.2 The plots show the reachable set computed for varying initial surge

velocities ranging from 1-5 m/s. The final plot shows the classification boundary with different initial conditions. . . 31 5.3 The plot shows the reachable sets for the marine vessel in the presence

of uncertainties. The uncertainty is present in the inertia matrix, where it can vary 25% of the Mass matrix. . . 31

vi

(9)

List of Figures vii

5.4 The plot shows the reachable set in the presence of disturbance. The disturbance in the direction of Surge, Sway and Yaw directions is assumed be within 1/10^ththe maximum actuator force. . . 32

(10)

(11)

List of Acronyms

DOF Degree of Freedom

GPS Global Positioning System

DP Dynamic Positioning

HJBI Hamilton Jacobi Bellman Isaacs PDE Partial Differential Equation 2PBVP 2 Point Boundary Value Problem CORA COntinous Reachability Analyzer

NED North-East-Down

CoG Center of Gravity

SVM Support Vector Machine

SQP Sequential Quadrartic Programming ODE Ordinary Differential Equations

ZOH Zero Order Hold

NLP Nonlinear Programming

ACADO Automatic Control And Dynamic Optimization PCA Principle Component Analysis

ix

(12)

(13)

Chapter 1 Introduction

1.1 Background

The once conservative maritime industry is now moving towards digitalization and autonomy. The advantages provided by increased autonomy in the maritime industry allows for lower fuel costs and also improves crew, passenger and cargo safety.

The change is largely brought about by the scrutiny of regulatory and environ- mental bodies for the shipping industry to be more efficient, modern and most importantly safe. In order to bring about autonomy, the first step is to provide situational awarenessto the captain, where he can perceive what is around the vessel and know about how his ship can maneuverer in the particular sea state. The next step in automating marine vessels would be to develop algorithms that act as decision support systems. Since the evolution of autonomy is a continual process, it is important to create algorithms that suggests the captain of what could be done.

The final goal is to develop algorithms that can make decisions on its own, but which can be overridden by the captains commands.

Whatever the level of autonomy in the marine vessels, the most important fea- ture that cannot be compromised upon is the guarantee of safety for the vessels.

The present work tries to solve the part of the puzzle, wherein the captain un- derstands his ship’s maneuverability better, thereby enabling him to take a safe course.

Many real world systems evolve according to complex nonlinear dynamics. The behavior of such systems is hard to predict and can be non-intuitive. Many of these systems are safety critical, which means verification plays an important role in checking the specification properties of the system. Verification of safety critical systems is challenging due to multiple reasons. The challenge stems from the fact that, the systems considered are generally nonlinear in nature, and also is high- dimensional. Further for a safety critical system, all the system behaviors must be accounted for, in the presence of disturbances and uncertainties in modeling.

Since there exists infinitely many possible trajectories that the system can evolve 1

(14)

2 CHAPTER 1. INTRODUCTION

with a given initial control and disturbance space, simulation based approaches are rendered useless. Thus formal verification methods are needed for guaranteeing performance and safety properties of systems.

Reachability analysis of marine vessels online, is a recent problem, where one tries to compute online the reachable set in the state space of the marine vessel in some amount of time in the future, given the input constraints due to the ship’s actuator system, and also given some bounds on the external disturbance. The external disturbance arises from various factors including, wind, wave and current forces, which cannot be measured accurately. The problem is nontrivial since the computation of the reachable set is difficult. The complexity for computation stems from:

• Complex nonlinear dynamics involving Coriolis force and nonlinear damping for the vessel.

• Maneuvering model is 3 DOF, which means 6 nonlinear ordinary differential equations are needed to characterize the behavior, which contributes to the curse of dimensionality.

1.2 Applications

The applications for reachability analysis for the marine vessels, extend the level of autonomy of vessels in general. The primary application of reachability for the marine vessel is for the captain to understand where the ship can maneuver. This offers a type of internal perception, for the captain to know about his ship’s ability to maneuver. Also, when the ship is in the open sea, there are quite long periods of time, the ship does not receive GPS coordinates, and the control systems requiring GPS data get affected. Reachable set could then be used as dead-reckoning for the controllers (or for manual control).

The next example extends to safety verification, as mentioned earlier. The reachable set entails all the possible state-space configurations, the dynamics can evolve in. This can be checked against an obstacle, to verify if the system is safe (Figure 1.1). Since, the algorithms developed in this thesis also provide the control signal along with the reachable set, a collision-avoidance algorithm can be formulated as a natural extension.

The area of reinforcement learning and the synthesis of learning based controllers, try to explore the state-space, to optimize a reward function to obtain a control signal. A learned control system may be guaranteed safe during deploy- ment, but it need not be safe while learning. Consider the example of a quadcopter learning to fly on a desired height. It needs to explore the state-space, and there might exist a search, where the quadcopter might reach an unsafe set of states (hitting the ground or the ceiling). So, reachability analysis can be used to ensure safety in this learning process.

(15)

1.3. OBJECTIVE 3

Figure 1.1: Verification of the vessel’s safety using the Forward reachable set.

As a final example, reachability analysis can be used to initiate parking of marine vessels. In the present day situation, the captain uses a ad-hoc method from experience, to understand when the ship needs to de-accelerate for parking. In many cases, this results in a not so smooth velocity profile for the ship. Backward Reachability analysis, answers the question of where should the captain initiate parking, if he needs to land in the parking spot at the present time.

1.3 Objective

The research objective is to compute online, the reachable set for the nonlinear ship maneuvering model, in the presence of disturbance in the form of waves, winds and currents and also taking into account uncertainties in the maneuvering model. The questions that are answered by the thesis are,

1. Given the azipod actuator constraints of the ship, where could the ship be?

(in a certain amount of time in the future).

2. Given the bounds for maximum allowable disturbance and uncertainty in the

(16)

4 CHAPTER 1. INTRODUCTION

vessel’s model, where could the ship be? (in a certain amount of time in the future.)

Since reachability for a marine vessel is a new problem, several methods were used to solve it. Each method reuires a trade-off between the generality of the dynamical system considered and the computation time. The most applicable so- lution, to our problem lays an emphasis on ’real-time computation’, which involves machine learning for the most approximate solution for the reachability problem.

1.4 Outline of the thesis

This thesis is organized in four chapters as follows. The first chapter motivates the reachability problem for marine vessels. Chapter 2 explains the state of the art in marine dynamics and how the disturbances & uncertainties are handled in the literature and also how reachable sets are computed. Chapter 3 explains nonlinear dynamic model used for maneuvering of the marine vessel. It also sum- marizes the preliminary results obtained using the Hamilton Jacobi based method and set-theoretic based methods. Chapter 4 talks about the major contribution of the thesis, which talks about the Classification based computation. Chapter 5 details the simulation results obtained for the marine vessel, in different experiments.

Finally, in Chapter 6, conclusions and future work are discussed.

(17)

Chapter 2 Literature Survey

Model based control for steering, dynamic positioning of ships has become state- of-the art since Linear Quadratic Regulators came to be used in the ship industry.

The exact model derived using the first principles, are given by a set of complicated differential equations describing 6 degrees-of-freedom (6 DOF) motion. For slow speed applications, such as Dynamic Positioning (DP), mooring and slow speed reference tracking, a simpler model can be considered, that is linear in the kinetic part. DP models and its simplifications are detailed in [4]. On the other hand, in high speed applications such as automatic course control, high speed position tracking, and path following, the nonlinear effects from the Coriolis and the damping terms become dominant. This model is developed in [5], where complete modeling, identification, and control design for maneuvering a ship along a desired path is presented. The identification is performed on a scaled down ship called CyberShip II. This thesis uses the model proposed in [5] for the computation of the reachable set.

Although Reachability analysis is a well-studied field in the control and the computer aided verification community, the problem is solved exactly only for a class of problems. Reachability analysis is studied extensively in computer software and in discrete state spaces, where properties can be formulated by the use of Linear Temporal Logic and other advanced formalisms that checks whether the runs of a transition system, or words of a finite automaton satisfy a set of properties [6].

This logic can specify the dangerous states the discrete state space needs to avoid, thereby verifying the safety properties of the system. While the reachability problem can be solved exactly and efficiently for discrete state spaces, the computation for continuous state spaces with more dimensions generally become increasingly difficult.

In [7], the Hamilton Jacobi method is developed, where it is shown that the solution can be obtained by computing the zero-level subset of the solution of Hamilton- Jacobi-Bellman-Isaacs Equations (HJBE). Unfortunately, while this methods can find accurate approximations to the reachable set for systems with complicated

5

(18)

6 CHAPTER 2. LITERATURE SURVEY

nonlinear dynamics, their computational cost scales exponentially with the system’s state space dimension 2.1. The major advantage of using HJ reachability stems from the fact that, the solutions are accurate without approximations, which means it can be used for guaranteed safety analysis.

Figure 2.1: Exponential increase in complexity as the number of dimensions in the state space increases [1].

Recently there have been several advances in HJ reachability theory and applications. In [8], a new technique is developed that decomposes the dynamics of nonlinear systems into subsystems which are coupled through common states, con- trols and disturbances, to tackle the curse of dimensionality. The reachable sets are computed for individual subsystems and the reachable set for the complete system is reconstructed without incurring additional approximation errors. Computations done using this techniques become degrees of magnitudes faster. In the case involving adversarial disturbances, their technique provides slighly conservative reachable sets. The authors in [9] try to approximate solutions of the HJ PDE by implement- ing and analyzing learning-based algorithms to approximate the solution of certain types of HJ PDEs using neural networks. The HJ reachability is applied to motion planning for a vehicle in [10], using an algorithm called FasTrack, which provides a safety controller for the vehicle along with a guaranteed tracking error bound.

The tracking bound is computed offline using the exact HJ reachability computation.The algorithm works by alternating between the target set being reach goal

(19)

7

set and the reach avoid set in the case of obstacles. In [11], a general framework is developed for ensuring safety while learning, implemented in quadcopters, which reduce the conservatism online as the system learns more information about the environment.

The intractability posed by the HJ reachability analysis calls for other approximate methods that rely on over approximate computation of the reachable sets, that check if the over approximate forward reachable set intersects a set of undesired states. This work uses the CORA toolbox [2] , which accepts continuous, nonlinear systems for computing the over approximate reachable set. It can also determine the reachable set for systems with parametric uncertainty. The sets are represented by zonotopes, and the error incurred due to linearization is used for enlarging the reachable set. The downside of this approach is that there is a trade-off between the accuracy of the reachable set, the nonlinearity in the system dynamics and the size of the admissible control input set. So, the over approximation of the reachable set becomes impractical, if both the nonlinearity and the admissible control set is large. Other tools such as SpaceEX [12], F low^∗ [13], can also can represent non-convex sets and are especially promising for reachability analysis of nonlinear models. The ellipsoidal toolbox [14] on the other hand approximates all the sets as a representation of ellipsoids, can be used to compute the reachable set of linear, continuous and hybrid systems.

The work done in [15], forms the basis for the main results derived in this thesis. In this work, supervised machine learning techniques are used to accurately predict cost-limited reachable sets of dynamical systems in real-time. A 2 Point Boundary Value Problem (2PBVP) is solved over the state space grid to create the dataset for the classification algorithms. The 2PBVP solver employs Cheby- shev Pseudospectral Methods and sequential quadratic programming. Two machine learning algorithms, Locally-weighed regression and Support Vector Machines, are used to demonstrate the method for query based classification of the reachable set.

The classical dubin’s car is used as a case study and the computation time required is 4 orders of magnitude less than the traditional exact approach. This technique is used, by the same authors to solve the real-time problem of kinodynamic motion planning for fixed-wing UAV navigating through a forest in [16]. The framework relies on a lookup table that stores precomputed optimal solutions to the 2PBVP.

(20)

(21)

Chapter 3 Problem Formulation

In this chapter, the details of the nonlinear maneuvering dynamical model for the vessel is developed. The model also incorporates the dynamics and the constraints posed by the azipod actuators of the vessel. Further, the reachability problem is posed and tackled using the Hamilton Jacobi Reachability method and by Set- theoretic computation schemes. The advantages and drawbacks of both the methods are evaluated.

3.1 Nonlinear Maneuvering Model of the Marine Vessel

When the ship maneuvers in the open sea, the ship can have motion in 6 degrees of freedom (DOFs). For the purposes of this this thesis, it is interesting to just consider the motion in the horizontal plane, characterized by Surge (Longitudinal forward motion), Sway (Sideways motion) and Yaw (rotation about the vertical axis) describing the heading of the vehicle. The earth-fixed reference frame {e}

is approximated by the North-East-Down (NED) frame. Let another coordinate frame {b} be attached to the center of gravity of the vessel. The discussion follows from the model posed by Fossen [5]

The kinematic relationship between the earth-fixed position vector η = [x, y, ψ]⁰ and the body-fixed velocity vector ν = [u, v, r]⁰ is given by,

˙η = R(ψ)ν (3.1)

where,

R(ψ) =





cos(ψ) −sin(ψ) 0 sin(ψ) cos(ψ) 0

0 0 1





Deriving the equations for the rigid body dynamics, using Newton’s second law, we arrive at,

MRB˙v + CRB(v)v = τ (3.2)

9

(22)

10 CHAPTER 3. PROBLEM FORMULATION

Figure 3.1: Motion in 6 degrees of freedom

where MRBis the rigid-body system inertia matrix, and CRBis the correspond- ing Coriolis matrix, and τ = [X, Y, N]⁰ is a generalized vector of Forces (X,Y) and Moment N. This is provided by the superposition of the actuator forces and moments. The external actuator forces are also represented in the body-fixed reference frame. Here, MRB and CRB take the form,

MRB





m 0 0

0 m 0 0 0 Iz



; CRB(v)





0 0 −mv

0 0 mu

mv −mu 0





The forces encountered by the ship’s hull due to its interaction with the water can be predominantly characterized by the hydrodynamic added mass, radiation- induced potential damping and restoring forces. Hydrodynamic potential theory programs can be used to compute the added mass and damping matrices by in-

(23)

3.2. THE REACHABILITY PROBLEM 11

tegrating the pressure of the fluid over the wetted surface of the hull or it can be identified during maneuvering experiments. For 3 DOF maneuvering, the restoring forces (hydrostatics due to the Archimedes principle) can be neglected. Here, the added mass (MA) and the corresponding Coriolis matrix (CA) take the form,

M_A=





−X_u_˙ 0 0

0 −Y_v_˙ 0 0 0 −N_r_˙



; CA(v) =





0 0 Y_v_˙v 0 0 −X_u_˙u

−Y_v_˙v X_u_˙u 0





The damping forces and moments can be represented by linear terms [Xu, Y_v, N_r].

For a constant surge speed, only these linear terms come into picture. But since we are seeking a globally valid model, a nonlinear representation called the second order modulus model where the coefficients are based on experimental data and curve fitting. Hence the total damping matrix can be collected into D(v) as,

D(v) =





−X_u− X_|u|u|u| 0 0

0 −Y_v− Y_|v|v|v| − Y_|r|v|r| −Y_|v|r|v| − Y_|r|r|r|

0 −N_|v|v|v| − N_|r|v|r| −Nr− N_|v|r|v| − N_|r|r|r|



 The final dynamics of the marine vessel can be written down as below, using superposition,

(MRB+ MA) ˙v + (CRB(v) + CA(v) + D(v))v = τ + w(t) (3.3) where w(t) represents a time-varying disturbance that is caused due wind, waves and ocean currents.

The actuators can be modeled using a first order system that models the time delay of the mechanical azipods to achieve a given force (X,Y) and moment (N) .

˙Ti= ui

Ts

− Ti

Ts

, where i = X, Y, N (3.4)

Hence the system is modeled using nine states x = (x, y, ψ, u, v, r, Tx, T_y, T_r) and three control inputs u = (ux, uy, ur). The concise notation for states and inputs being, x = (η, ν, τ) and u respectively.

3.2 The Reachability Problem

In the reachability problem, we try to determine the set of states the system can reach, under the influence of a set of control inputs & disturbances and a set uncertain parameters. This work concentrates on three different approaches for computing the reachable set for the vessel.

• Hamilton Jacobi Reachability Analysis

(24)

• Set-Theoretic Reachable set computation.

• Classification based Reachability Algorithm

The next sections, in this chapter, detail the first two methods, its application to ship’s reachability analysis and its drawbacks. The next chapter details the main results of the thesis for computing the reachable set for the marine vessel, using the third method, where supervised machine learning is used to compute the reachable set and how the drawbacks posed by both Hamilton Jacobi analysis and the set- theoretic approaches are overcome.

3.2.1 Hamilton Jacobi Reachability Analysis

Hamilton Jacobi reachability analysis falls under the umbrella of optimal control problems and differential games, which are theoretical tools that can be used for verification of safety-critical systems. The reachable set of a dynamical system can be computed with respect to a target set, which describes the set of final conditions under consideration. The final goal could be a desirable set of states (goal set) or an unsafe set of states (avoid set). The goal of Hamilton Jacobi (HJ) reachability is to compute different flavours of reachability namely, Backwards Reachable set and the Forward reachable set. Let us setup the mathematical theory behind HJ reachability according to [7].

Figure 3.2: Illustration of the level set method, which encodes the target set and the reachable set using value functions.

(25)

Let x ∈ Rⁿ be the system state, which evolves according to,

˙x = f(x(t), u(t), w(t)), u(t) ∈ U, w(t) ∈ W (3.5) where, u(t) is the control input and w(t) is the disturbance acting on the system.

The sets W and U represent the control and state constraints respectively. Intu- itively, a Backward Reachable set represents the set of states x ∈ Rⁿfrom which the system can be driven into some set <0⊆ Rⁿat the end of time horizon t. Here, the control input tries to steer the system towards the goal, whereas the disturbance tries to steer away from the goal. Hence, we want to compute the set,

<(t) = {x : ∃w(·) ∈ W, ∀u(·) ∈ U, x(0) ∈ <0} (3.6) Formally, let Jt(x, u(·), w(·)), be the cost accumulated, after a finite time horizon.

If we consider Jt(·) as the signed distance between the system state and the target region at the terminal state of the system, we can determine if the trajectory reached the target if the signed distance is negative. The target set, <0can be encoded as the zero sub level set of a bounded and Lipschitz continuous function v(x), such that,

<0= {x ∈ Rⁿ|v(x) ≤ 0} (3.7) The reachable set can be encoded, using the value function as,

<(t) = {x ∈ Rⁿ|V(t, x) ≤ 0} (3.8) If we define the cost as,

Jt(x, u(·), w(·)) = v(x(0)) (3.9) then the system reaches the target set under the control u(t) and disturbance w(t) only if, Jt(x, u(·), w(·)) ≤ 0. The idea of the level-set method is illustrated in the Figure 3.2.

Since the control input wants to drive the system to the target, it wants to minimize the cost in 3.9, whereas disturbance wants to maximize ( or steer away from the target set) the cost. Hence, there is a dynamic game played between the control input and the disturbance input, resulting in,

maxw min

u J_t(x, u(·), w(·)) (3.10)

Solving 3.10, using the principle of dynamic programming, we land up in the Hamilton-Jacobi-Isaacs equation,

∂V(t, x)

∂t + H(t, x, λ) = 0, V(0, x) = v(x) (3.11) H(t, x, λ) = max

w min

u λ.f(x, u, w) (3.12)

(26)

The interpretation of <(t) is that, if x(t) ∈ <(t), then the control input can drive the system to the target set at time 0, irrespective of the disturbance. Here, if the target state is desirable, then we are interested in verifying if there exists a control input that can steer the system to the target despite the worst case disturbance. But if the target set <0is an obstacle (avoid set), we need to minimize over disturbance and maximize over the control input.

3.2.2 Implementation

The implementation of Hamilton Jacobi Reachability analysis in Level-Set toolbox[17]

in Matlab, using the ship’s dynamics, was carried out. Since the computation time grows exponentially with the number of states [7], A simpler system, taking into account only the surge and Sway velocity states, was implemented in the Level-Set toolbox. Referring to the surge and sway equations (equation 3.3), it can be seen that the control (u(t)) and the disturbance inputs (w(t)) enter the system affinely.

The optimal control input and the optimal disturbance input is required to compute the Hamiltonian required to solve the PDE (equation 3.11). If the control inputs are constrained (as in our case), such that u(t) ∈ [umin, umax], w(t) ∈ [wmin, wmax] this results in a bang-bang control. The inputs to the toolbox are the functions representing,

•Dynamics equations in ODE

˙x = f(x(t), u(t), w(t)) = F (x(t)) + u(t) + w(t) (3.13)

•Optimal control input

u^∗(t) = argmaxu{λ.f(x(t), u(t), w(t))} = argmaxu{λ.u}=

(u_max, if λ > 0 umin, if λ ≤ 0 (3.14)

•Optimal disturbance input

w^∗(t) = argminw{λ.f(x, u, w)} = argminw{λ.w}=

(wmin, if λ > 0 wmax, if λ ≤ 0

(3.15) Figure 3.3, shows the forward reachable set, where the target set is a goal set, hence

(27)

the control tries to maximize the Hamiltonian and the disturbance tries to minimize the Hamiltonian.

The reachable set is computed for 1 s forward in time, but the computation for

Figure 3.3: Forward Reachable set for the velocity states in the surge and the sway direction.

solving the HJ equation takes 77s. When including the third state, the computation takes minutes, and if all the six states, including the kinematic and dynamic relationships, the program becomes intractable, and the Level-set toolbox throws an error statement that, "Dimension > 5 is dangerously large". Hence, although this method, for computing the reachable set provides lots of advantages, such as, providing theoretical guarantees, ability to handle sets of arbitrary shape and providing the optimal control inputs for reaching the desired state (or avoiding, based on the application), Hamilton Jacobi reachability suffers from the curse of dimensionality, and proves to be incompatible for the application at hand.

3.2.3 Over-approximate Set Theoretic Reachable set computation

The exact reachable set computation, as seen in the previous section, is time con- suming and proves intractable as the number of dimensions in the state space increases. Hence, an over-approximation of the reachable set can be employed, where the accuracy of the reachable set is compromised for faster computation time. Over-approximation of the reachable set is useful, in the context of safety verification. Clearly, if the over-approximated set of reachable states does not in- tersect the unsafe set, it means the original dynamical system is safe as well (Figure 3.4.

(28)

Figure 3.4: Illustration of the overaproximation of reachable set used for safety verification [2].

This thesis uses the Continuous Reachability Analyzer-CORA [2], which is a formal verification tool of cyber-physical systems. All the sets (in Rⁿ) are repre- sented, by a Zonotope representation, which are parametrized by a center c ∈ Rⁿ and generators gⁱ∈ Rⁿ, defined as,

Figure 3.5: Illustration of the steps involved in approximate reachable set calculation [2].

Z= {c +

p

X

i=1

βigⁱ|βi∈[−1, 1], c ∈ Rⁿ, gⁱ ∈ Rⁿ} (3.16) Zonotopes are a compact way of representing sets in high dimensions and, opera- tions required for reachability analysis, such as, linear maps, M ⊗Z := {Mz|z ∈ Z}

and Minkowski addition Za⊕ Zb := {za+ zb|za ∈ Za, z_b ∈ Zb} can be computed efficiently [2]. For a linear system, characterized by,

˙x = Ax + Bu x(0) ∈ X0⊂ Rⁿ, u(t) ∈ U ⊂ R^m (3.17) The idea is to first compute the the set of all solutions, from the initial set of states X0, which involves the computation of the matrix exponential e^At but assuming no uncertainty in the control input. The next step is compute a convex hull which encompasses the solution and the initial state. The final step is to enlarge the computed convex hull to account for the uncertain input (Figure 3.5).

(29)

3.3. SUMMARY 17

Since, the marine vessel is a nonlinear system, the above algorithm can be applied only if the nonlinear system is linearized. In CORA, the system is linearized in the center of the present state zonotope and the control set zonotope. The linearization error is accounted for, by including it in as an uncertain input, thereby enlarging the approximate reachable set. The Figure 3.6, shows the reachable set computed for the marine vessel using the kinematic and dynamical relationships (equation 3.1, 3.3).

Figure 3.6: Reachable set for the marine vessel with 3 DOF maneuvering model, for 3 min into the future, admissible control being 1/10^ththe original value.

The assumed constraints on the control inputs is one-tenth the original value for the real vessel. If the admissible control set is enlarged than this value, the uncertainty in the system explodes, thereby rendering an impractical over-approximation.

Since the nonlinearity in the dynamics of the marine vessel is huge, the admissible control input set cannot be high, because both contribute to the error in approximation of the reachable set. The computation of the forward reachable set for 3 min into the future, requires 5 s of computation time. In the Figure 3.6 , it can be seen that the solution obtained by this method is clearly the over-approximation of the true trajectories (black lines representing random trajectories).

3.3 Summary

In this chapter, the nonlinear maneuvering model of the ship is introduced, along with the actuator dynamics and constraints. The nonlinear model is used to compute the reachable set of the marine vessel, given the actuator limits and the limits

(30)

on the maximal allowable disturbance. The first proposed method using the theory of Hamilton Jacobi reachability analysis provides exact calculation of the reachable set, by solving the HJI PDE, using the level set method. But on the downside, the computation is O(eⁿ), where n is the number of states for the system. Hence, it proves intractable for the ship dynamics, which has as many as 6 states (without considering the actuator dynamics). In the other method involving, set theoretic computations, we try to linearize the nonlinear model and to enlarge the evolved reachable set with this linearization error. This method results in impractical reachable sets, as the linearization error grows, or if the admissible control range is large.

The results obtained in this chapter, using both the Hamilton Jacobi analysis, and the set theory based computation methods, motivated the author to use a supervised Machine learning based tequenique to best approximate the reachable set. This forms the content of the next chapter.

(31)

Chapter 4 Classification based Reachability Analysis

This chapter focuses on the implementation of the classification based algorithm, used for computing the reachable set for the marine vessel. The first step in the

Figure 4.1: Simplified 2D illustration of a cost limited reachable set for a given threshold cost of J0

algorithm is to sample the state space randomly and combine two points which act as the start point and the goal point. The next step is to compute a time optimal control strategy which minimizes the final time, where the trajectory is constrained to start and finish at particular points in state-space. If the time required to reach the goal point, is less than the threshold final time as required by the user, we

19

(32)

20 CHAPTER 4. CLASSIFICATION BASED REACHABILITY ANALYSIS

label the start-goal pair to be reachable, with the threshold final time (Figure 5.4).

This is iterated for different start-goal pairs and each pair is labeled to be either reachable or non-reachable. Support Vector Machines is used for separating the boundary between the reachable and non-reachable points. Given a new query start-goal pair, the trained Support Vector Machine can classify the query point as either reachable or non-reachable. This work is based on the theory formulated in [15] and builds on it for the computation of reachable sets in the presence of uncertainties and disturbances.

The organization of this chapter is as follows. Section 4.1 develops the optimal control formulation. In Section 4.1.1, we solve the optimal control problem using Multiple Shooting method and Sequential Quadratic Programming (SQP).In Sec- tion 4.2 it is shown how the developed dataset can be used to learn the nonlinear boundary separating the reachable and the non-reachable points in space. Later, in Section 4.3.1 and Section 4.3 talks about creation of the dataset and how it is implemented with the ACADO Toolbox [18] respectively Finally, in Section 4.3.2 an extension is presented to include uncertainties and disturbances into the solution framework.

4.1 Optimal Control Problem

Consider the nonlinear model of the ship, including the actuator dynamics (equation 3.3,3.4) represented by ˙x = f(x(t), u(t)), where x ∈ which evolves over a time horizon [t0, t_f].

Let us define a cost,

J= φ(x(t), tf)

| {z }

Mayer term

+Z tf

t0

L(x(τ), u(τ), τ)dτ

| {z }

Lagrange Term

(4.1)

and the corresponding cost-limitted reachable set, for the arbitrary start state xa as [19]

R(xa, U, J0) = {xb∈ Xk∃u ∈ U & ∃t⁰ ∈[t0, tf] (4.2) s.t. x(t⁰) = xb & J^∗≤ J0} (4.3) where x(t) and u(t) represents the state vector and control vector respectively, X represents the state constraints and U represents the admissible control inputs. J^∗ is the optimal cost for going from start (xa) to the goal (xb) state. Equation 4.2 represents the set of all goal states, which obey the state constraints, such that the optimal control responsible for arriving at the goal state obeys control constraints and that the cost for arriving at the goal state is less than the threshold cost J0. Since, we are interested in the minimum time optimal control for the vessel, the cost J = tf needs to be minimized. Specifically the time optimal solution J^∗ is

(33)

4.1. OPTIMAL CONTROL PROBLEM 21

given by:

J^∗=minimize

u

Z tf

t0

dt

subject to u(t) ∈ U (actuator Constraints)

x(t) ∈ X (state Constraints)

˙x(t) = f(x, u) (Ship Dynamics + Actuator Dynamics) x(0) = xa, x(tf) = xb

(4.4) Equation 4.18 is referred to as 2 Point Boundary Value Problem (2PBVP), trajec- tory optimization problem or the steering problem.

4.1.1 Optimal Control Solution

The first step in solving the optimal control problem is transcription. Since the optimization presented in 4.18 is continuous, it has infinitely many decision variables. Thus the problem 4.18 needs to be discretized. Multiple shooting method is used for discretization. The discretization leads to a Nonlinear program which is solved using a Sequential Programming Method. Following is the details of the multiple shooting algorithm and Sequential Quadratic Programming for solving the proposed optimal control problem.

Figure 4.2: Illustration of the multiple shooting algorithm [3]

4.1.1.1 Multiple Shooting

Multiple shooting works by breaking up a trajectory into some number of segments, and guessing the solution in these grid points Figure 4.2. The solution approach is analogous to shooting a cannon at a target. The algorithm first simulates the trajectory of the cannon and tries to minimize the defect between the shot, final

(34)

cannon state and the target state. The simulation is usually broken down into multiple segments (Multiple Shooting) and it is stringed together by the optimization algorithm.

The continuous control input is discretized using a polynomial in the general case or using Zero Order Hold (ZOH) in the simple case in the fixed time grid 0 = t0 < t₁ < ... < tN = T . The discrete control represented by ZOH, say u(t, q) = qk for t = [tk, tk+1]. The ODE, present in the optimal control problem is simulated separately on each iteration interval [tk, tk+1], using the ship dynamics, starting with artificial initial values si:

˙xk(t, sk, qk) = f(xk(t, sk, qk, qk)), t ∈ [tk, tk+ 1] (4.5)

x(tk, s_k, q_k) = sk (4.6) Similarly, we numerically compute the integrals,

l_k(sk, q_k) =Z tk+1

t_k

L(xk(tk, s_k, q_k), qk)dt =Z tk+1

t_k

dt (4.7)

The final issue of stringing together of the trajectories is left to the NLP solver as constraints. Also, we choose the time grid on which the inequality path constraints are checked (commonly the same grid as the one used for control discretization).

So we arrive at the final NLP formulation,

minimize_s,q φ(sN) +

N −1

X

i=0

l_i(si, q_i)

subject to x0= s0 (x0= xa) (initial value)

x_i(ti+1, s_i, q_i) = si+1≤0 i = 0, ...N − 1 (Stitching the trajectories) h(si, qi) ≤ 0 i = 0, ...N − 1 (path constraints)

r(sN) = 0 (xN = xb) (terminal constraints) (4.8) 4.1.1.2 Sequential Quadratic Programming

The nonlinear program is solved by linearizing the cost function and the constraints.

This yields a quadratic programming subproblem. The SQP problem requires the cost function, and the constraints to be twice continously differentiable, since a Hessian needs to be computed at each time step. The nonlinear program (Equation 4.8) obtained through multiple shooting can be summarized as,

minimize

x h(x)

subject to g(x) ≤ 0 r(x) = 0

(4.9)

(35)

4.2. SUPPORT VECTOR MACHINE FOR BOUNDARY CLASSIFICATION23

The Lagrangian of the problem is defined as,

L(x, λ, σ) = h(x) + λ^Tg(x) + σ^Tr(x) (4.10) where, λ and σ are the lagrange multipliers. A guess point xk is chosen and a search direction dk is chosen, which is the solution of the quadratic programming subproblem,

minimize

d h(xk) + ∆h(xk)^Td+1

2d^T∆²xxL(xk, λ_k, σ_k)d (4.11)

subject to g(xk) + ∆g(xk)^Td ≤0 (4.12)

r(xk) + ∆r(xk)^Td= 0 (4.13)

4.2 Support Vector Machine for boundary classification

The idea behind this approach is to approximate the reachability boundary using a nonlinear binary supervised classification algorithm [20]. It tries to separate the training (start, goal) pairs as being reachable or non-reachable sets. The trained nonlinear classifier can then be used to predict whether a new pair of (start, goal) pair is reachable or not. For the application of computing the reachable set for a marine vessel, in this thesis, a binary Support Vector Machine algorithm with a nonlinear kernel function is employed. SVM seeks to design a function h(w^Tx+ b) affinely dependent on the training dataset ˜x with the weights w and a bias term b, that correctly maps a vector ˜x to a label y ∈ {−1, +1}, which is the separating hyperplane between the two classes. In our example, -1 means that the given (start, goal) pair is not reachable within the defined cost threshold, whereas +1 means the pair is reachable. The parameters are designed to minimize the misclassification that occur in the training set, without over fitting the data. This optimization problem can be formulated as: finding b and w, that minimize kwk, such that for all the data points ( ˜xi, yi),

y_ih(xi) ≥ 1 (4.14)

The support vectors are the xi on the boundary, those for which yih(xi) = 1.

The resulting problem is a quadratic programing problem, and the optimal solution ( ˆw, ˆb) enables the classification of a new test data point z as,

class(z) = sign(˜zˆw+ ˆb) (4.15)

It is often computationally simpler to solve the dual of the above problem by taking the positive Lagrange multipliers αj multiplied to each constraint and sub-

(36)

Figure 4.3: Illustration of separating hyperplane & support vectors

tracting from the objective function to get,

maximize_α

m

X

i=1

α_i−1 2

m

X

i,j=1

yⁱy^jα_iα_jh ˜xⁱ,˜x^ji subject to αi≥0, i = 1, ..., m

m

X

i=1

αiyⁱ= 0

(4.16)

In the above case, it is assumed that the data is perfectly separable, but in reality l¹-regularization is performed which introduces slack variables to 4.16. Also since we are dealing with a nonlinear boundary classification problem, there is no simple hyperplane separating the two types of data. A kernel function creates linearly separable data in the higher dimension. Hence a kernel function (Gaussian in this case), k(˜x, ˜z) = e^k˜^x−˜^zk² replaces the dot product in the above set of equations.

4.3 Implementation in ACADO

The previous sections introduced the general theory behind the solution of the optimal control problem. In this section, we will adapt the theory to solve the time- optimal control of the ship, including the actuator dynamics. Since the SQP method introduced in Section 4.1.1.2, requires the constraints to be twice differentiable, and the ship dynamics (Equation ), pertaining to the velocity states have absolute value functions, an approximation needs to be carried out.

(37)

4.3. IMPLEMENTATION IN ACADO 25

Approximation 1: The cubic polynomial fit for the ν|ν| term in the velocity states, i.e.

ν|ν|= a0ν³+ b0ν²+ c0ν+ d ν = [u, v, r]⁰ (4.17)

Approximation 2: The cross coupling terms, |v|r, |r|v in the differential equations pertaining to sway and yaw rate dynamics are assumed to be small.

The cubic approximations are carried out in the state space domains of the individual velocity states, providing a really good fit to the actual dynamics. The above approximations work quite well in converting an ill-posed optimization problem to a problem that can be solved robustly. Moreover, the cross-coupling terms in the sway and yaw rate dynamics can be assumed to be zero while considering the case of maneuvering, as neither sway nor yaw rate is high during maneuvering in the open sea.

In this thesis, the ACADO Toolkit [18], which is a software environment and algorithm collection for automatic control and dynamic optimization has been used to solve the optimal control problem:

J^∗=minimize

u T

subject to u(t) ∈ U (actuator Constraints)

x(t) ∈ X (state Constraints)

˙x(t) = f(x, u) (Ship Dynamics + Actuator Dynamics) x(0) = xa, x(tf) = xb

(4.18) Here, the differential equation constraint corresponds to the a kinematics and the dynamics of the ship (Equation [3.3]), as well as the dynamics of the actua- tor(Equation [3.4]) (modeled by the first order system in each DOF). U corresponds to how the ship’s azipods are constrained by the maximum and minimum thrust that can be produced in all the 3 DOF. The states corresponding to heading and the velocity states are also constrained by a maximum value, when the ship is maneuvering in the open sea.

Problem 4.12 can be cast in ACADO [18], which tries to compute the time- optimal control policy for the ship starting from xa to reach xb. ACADO code is written so that Problem 4.18 is transcribed by the Multiple Shooting technique. The control vector is discretized using 40 discretization points. The maximum final time threshold is 200 s, 0 ≤ T ≤ 200. So the discretization interval is 5 s. For the above posed problem with box input constraints, the exact Hessian can be computed, which can give results with higher accuracy, albeit being computationally slow.

Since we are dealing with a nonlinear optimization problem, the sufficient conditions of optimality can only guarantee the convergence to a local minima. The exact Hessian works in our favour, where the solution almost always converges to the global minima.

(38)

4.3.1 Dataset Creation

The dataset that needs to be trained by the supervised binary classifier requires two labels of either being reachable or not. The dataset is created by changing different values in the optimization problem and solving it to gather the information of it being reachable or not.

The important details of how the initial and the final values of the different states are constrained, according to the condition of maneuvering in the open sea, are outlined as follows. The initial values of states x, y, ψ can without loss of generality be chosen as zero, since the acquired reachable set can be translated to any x, y and be rotated with the current heading ψ to obtain the reachable set at any x, y, ψ. Among the velocity states (u, v, r), during maneuvering, it can be assumed for the sake of minimizing the complexity, only the surge velocity changes (from zero to maximum allowed surge speed u = umax, with linear interpolation in between), and the sway velocity and yaw rate remains close to zero, since, during maneuvering, the ship is not expected to turn unexpectedly. The initial state for the thruster in the surge direction is constrained by the steady-state thruster value, which has been responsible for producing the initial surge velocity (Line 7 in Algorithm 1). The final states of (x, y) take the form of a sampled grid in space in front of the marine vessel. All the other states are not constrained, and are left to take any value.

To summarize, the problem 4.18 needs to be solved multiple times looped over different initial surge velocities (interpolated between u = 0 to u = umax), and also over the grid of final spacial states of (x, y). Finally in the online phase, the present state information of the ship is received (Line 17 in Algorithm 1), and the SVM predicts the reachable set boundary for the given query surge velocity. In order to plot the boundary, various final state points are plotted. This boundary represents the reachable set, if the ship was at poisition η = [0, 0, 0]⁰. So the position states needs to be translated and rotated with the current position information (Line 19,20 in Algorithm 1). Finally, the rotated and translated reachable set is plotted.

(39)

4.3. IMPLEMENTATION IN ACADO 27

Algorithm 1:Classification based Reachable set Computation

1 OFFLINE

2 input: u_max, T_threshold

3 Interpolate: uvector = [0, umax]

4 Interpolate: xf inal = 2-D Mesh in (x,y);

5 while x_b ∈ x_{f inal} && uinit∈ u_vector do

6 T hrust0 = Solve steady state surge equation with [ ˙u, v, r] = 0 & u = uinit;

7 T_opt= Solve Problem 4.18 with

xa= [0, 0, 0, uinit,0, 0, T hrust0,0, 0] & xb(1 : 2) = xf inal;

8 if Topt≤ T_threshold then

9 label=1;

10 else

11 label=-1;

12 end

13 dataset= [xf inal; uinit; label];

14 end

15 SV MM odel= Train(dataset);

16 ONLINE

17 x_now = xi, y_now= yi, ψ_now= ψi, u_now = ui;

18 [xreach, yreach] = predict boundary(SV MM odel(uinit));

19 x_solution = xnow+ xreachcos(ψnow) + yreachsin(ψnow);

20 y_solution= ynow− x_reachsin(ψnow) + yreachcos(ψnow);

21 Plot(xsolution, ysolution);

4.3.2 Including Uncertainties and Disturbances

The uncertainties in the system dynamics and external disturbances due to the sea condition alter the reachable set R. The idea of including uncertainties is as follows. The inertia matrix could be uncertain, due to the loading conditions on the vessel. Hence the inertia matrix can vary from [M − ∆M, M+ ∆M], where ∆M

corresponds to the uncertainty introduced due to the different loading conditions of the marine vessel. The algorithm presented in Section 4.3 can be carried out with varying inertia matrix (linearly interpolated from M − ∆M to M + ∆M). Let it produce m such varying Inertia matrices. We obtain m such reachable sets <1...<m. If the target set is a goal set (defined in ...) the overall reachable set corresponds to <uncertain= <1∩ <2.. ∩ <m. The overall reachable set is an intersection of the individual sets, since, for all the possible uncertainties in the marine vessel, there exists a control that can steer the ship to the goal set <uncertain. Also, for the computation of the avoid set or the unsafe set (defined in ...), the reachable set

(40)

becomes <uncertain= <1∪ <2.. ∪ <_m. The overall reachable set is the union of the individual sets, since, despite all the possible control actions of the marine vessel, there exists an uncertainty that can drive the vessel to the unsafe set <_uncertain.

A similar approach can be adopted to include disturbance, wherein the disturbance can act in all the 3 DOF in either the positive direction or the negative direction.

−d^max_u ≤ wu(t) ≤ d^maxu (4.19)

−d^max_v ≤ wv(t) ≤ d^maxv (4.20)

−d^max_r ≤ wr(t) ≤ d^maxr (4.21) where, w(t) = [wu(t), wv(t), wr(t)]⁰ represents the disturbance vector and the maximum values are represented by dmax = [d^maxu , d^max_v , d^max_r ]⁰. Since, the disturbance shrinks or enlarges the reachable set, either the maximum value or the minimum value (maximum in the negative direction) produces the most deviation from the normal reachable set(< is monotonous with respect to disturbance w(t)).

Hence there are eight possible combinations of possible disturbance values (either d_maxor −dmax). Similar to the approach used for uncertainties, the goal reachable set is <disturbance= <1∩ <₂.. ∩ <₈and the avoid set is <disturbance= <1∪ <₂.. ∪ <₈ .

4.4 Summary

In this chapter, the reachable set is computed in real time, in three major steps.

The first phase is the creation of dataset, wherein a time-optimal control problem is solved using various features such as varying initial surge velocity, variation in the inertia matrix, and varying disturbances in the 3 DOF, and if the optimal time is less than the threshold time required, it is labeled as +1 and if not, it is labeled as -1. The required reachable set (goal or avoid set) in the presence of disturbances and uncertainties can be calculated either using the intersection or the union of the different sets. The next phase of the algorithm is learning a nonlinear binary SVM classifier using the created training dataset to separate the boundary of points in state-space, where the ship can reach in a certain threshold time. The final stage is the prediction of the reachable set, using the state information and as the input to the trained SVM model, in real time.

(41)

Chapter 5 Simulation Results

In this chapter, the results for the classification based reachability analysis has been detailed. The model used for the simulations is briefed in Chapter 3. The method used for the computations is detailed in Chapter 4.

5.1 Maneuvering reachable set for different initial conditions

The plot 5.1 shows the reachable set for the marine vessel computed using the classification based method. The state space is sampled, with 50mX50m boxes.

The initial surge velocity is taken to be 1 m/s. The computation assumes, that the reachable set is computed for 200 s in the future. The red points are those with the label of feasibility for the optimal control problem, while the blue points are not feasible within the given time frame.

The plot 5.2 shows how the reachable set is different for different initial velocities. This can be used to compute the reachable set online, while the ship is moving, where the surge velocity is fed back from the sensors. The final plot in Figure 5.2 shows the classification boundary with different initial conditions.

29

Real Time Reachability Analysis for Marine Vessels

Real Time Reachability

Analysis for Marine Vessels

SUDAKSHIN GANESAN

Acknowledgements

Contents

List of Figures

List of Acronyms

Chapter 1

Introduction

1.1 Background

1.2 Applications

1.3 Objective

1.4 Outline of the thesis

Chapter 2

Literature Survey

Chapter 3

Problem Formulation

3.1 Nonlinear Maneuvering Model of the Marine Vessel

3.2 The Reachability Problem

3.3 Summary

Chapter 4

Classification based Reachability Analysis

4.1 Optimal Control Problem

4.2 Support Vector Machine for boundary classification

4.3 Implementation in ACADO

4.4 Summary

Chapter 5

Simulation Results

5.1 Maneuvering reachable set for different initial conditions