Dynamics and Simulation of a Snakeboard‟s Motion

Dynamics and Simulation of a Snakeboard‟s Motion Jonathan Jamieson

Department of Mechanical and Aerospace Engineering

Supervisor – Dr James Biggs

MEng Degree in Mechanical Engineering

April 2012

i

Acknowledgments I would like to thank Dr James Biggs for agreeing to supervise this project and pointing me in the right direction.

ii

Abstract The aim of this paper was to investigate and understand the non-intuitive dynamics of a snakeboard with a long-term vision of developing controls for similar systems. A numerical simulation of the four, coupled, non-linear differential equations that represent the motion of the snakeboard was created using the MATLAB coding environment. Sinusoidal, periodic controls were defined and by varying the ratio of angular frequencies, the following gaits could be analysed: resonant, rotation and parking. Through experimentation with controlled parameters, it was found that the direction of travel for the resonant and parking gait depended on the amplitude relationship between the rotor and the platforms. The direction of travel for the rotation gait depended solely upon the amplitude of the rotor. The variation in behaviour caused by changing the amplitude and absolute angular frequencies was also investigated. Larger angular frequencies resulted in greater displacements but the change in behaviour using different amplitudes was more complicated. A bifurcation study was also performed. An open loop control method was developed so changes to controlled parameters could be made in real time and their effects observed. To provide more control, a physical input device was created to interface with a computer running the simulation. A numerical motion planning algorithm was developed to find the control parameters which would displace the snakeboard from rest, for a specified displacement in the 𝑥-direction in a given duration. It was found that precise angular frequencies were needed for accurate displacement. Consequently, the iterative process must refine the parameters being changed for a reasonable computational time. Finally, an algorithm that ensured stability in the resonance gait was developed. The behaviour matched closely matched the author‟s experience of riding snakeboard.

iii

List of Contents Acknowledgments ...........................................................................................ii Abstract .......................................................................................................... iii List of Contents ..............................................................................................iv List of Notations ........................................................................................... viii 1.

Introduction ............................................................................................. 1 1.1. Nonholomonic Systems .................................................................... 1 1.2. Background to Snakeboarding .......................................................... 3

2.

Literary Review ....................................................................................... 5

3.

Mathematical Model ................................................................................ 6 3.1. Nonholomonic Constraints ................................................................ 6 3.2. Controlled Parameters ...................................................................... 7 3.3. Physical Parameters ......................................................................... 7 3.4. ODE45 Differential Solver ................................................................. 7 3.5. Dynamical Equations of Motion ........................................................ 9 3.6. Validation ........................................................................................ 10 3.6.1.

Validation Test 1 ................................................................... 10

3.6.2.

Validation Test 2 ................................................................... 11

3.7. Limitations....................................................................................... 12 3.8. Plotting and Visualisation ................................................................ 13 3.8.1.

Static Overlay........................................................................ 13

3.8.2.

Animated 2D ......................................................................... 14

3.8.3.

Animated 3D ......................................................................... 14

3.9. MATLAB Code ................................................................................ 14 4.

Resonance Gait .................................................................................... 15 4.1. Direction of Travel ........................................................................... 16 4.2. Amplitude ........................................................................................ 17

iv

4.3. Angular Frequency ......................................................................... 19 4.4. Phase shift ...................................................................................... 20 4.5. Instability ......................................................................................... 21 5.

Rotation Gait ......................................................................................... 23 5.1. Direction.......................................................................................... 24 5.2. Rate of rotation ............................................................................... 25

6.

5.2.1.

Angular frequency ................................................................. 25

5.2.2.

Amplitude of platforms and rotors ......................................... 26

5.2.3.

Initial Velocity ........................................................................ 27

Parking Gait .......................................................................................... 28 6.1. Direction.......................................................................................... 28 6.2. Rate of change in the 𝒚-direction .................................................... 30 6.2.1.

Angular Frequency ................................................................ 30

6.2.2.

Amplitude of platforms and rotors ......................................... 31

6.3. Initial Velocity .................................................................................. 32 7.

Direct Simulation Controller .................................................................. 33 7.1. Limitations....................................................................................... 34 7.2. Construction.................................................................................... 35 7.2.1.

Components.......................................................................... 35

7.2.2.

Microcontroller ...................................................................... 35

7.2.3.

Final Design .......................................................................... 35

7.2.4.

Improvements ....................................................................... 36

7.3. Results ............................................................................................ 36

8.

7.3.1.

Validation Test ...................................................................... 37

7.3.2.

Resonance imitation test ....................................................... 38

7.3.3.

Rotation Imitation Test .......................................................... 40

Open loop periodic controls .................................................................. 42

v

9.

Bifurcation study ................................................................................... 44 9.1. Phase shift ...................................................................................... 44 9.1.1.

Rotor phase shift ................................................................... 44

9.1.2.

Platform phase shift .............................................................. 45

9.2. Angular Frequency ......................................................................... 46

10.

9.2.1.

Rotor angular frequency........................................................ 46

9.2.2.

Platform angular frequency ................................................... 47

Motion Planning ................................................................................. 49

10.1. Speeding up and slowing down ...................................................... 49 10.1.1.

Test 1 .................................................................................... 49

10.1.2.

Test 2 .................................................................................... 50

10.2. Parameter to be changed ............................................................... 51 10.2.1.

Rotor Amplitude Test ............................................................ 51

10.2.2.

Platform Amplitudes Test ...................................................... 52

10.2.3.

Velocity Test ......................................................................... 53

10.2.4.

Conclusions .......................................................................... 54

10.3. Algorithms ....................................................................................... 55 10.3.1.

Version 1 ............................................................................... 56

10.3.2.

Version 2 ............................................................................... 58

10.3.3.

Version 3 ............................................................................... 59

10.4. Improvements ................................................................................. 61 11.

Stable Resonance Gait ...................................................................... 62

12.

Conclusions ....................................................................................... 66

12.1. Summary ........................................................................................ 66 12.2. Evaluation ....................................................................................... 67 12.3. Further Work ................................................................................... 67 13.

References ......................................................................................... 68

vi

Appendix I

– Timestep .............................................................................. 1

Appendix II

– MATLAB Code ..................................................................... 3

Appendix III

– Microcontroller Code .......................................................... 26

vii

List of Notations 𝐽𝑏

–

Moment of inertia of crossbar (kg·m²)

𝐽𝑝

–

Moment of inertia of platforms (kg·m²)

𝐽𝑟

–

Moment of inertia of rotor (kg·m²)

𝑎𝑏

–

Amplitude of rotation for the back platform (radians)

𝑎𝑓

–

Amplitude of rotation for the front platform (radians)

𝑎𝑟

–

Amplitude of rotation for the rotor (radians)

𝛽𝑏

–

Phase change of the front platform (radians)

𝛽𝑓

–

Phase change of the front platform (radians)

𝛽𝑟

–

Phase change of the rotor (radians)

𝜑𝑏

–

The back platform‟s angle of rotation with respect to the centre line of the snakeboard (radians)

𝜑𝑓

–

The front platform‟s angle of rotation with respect to the centre line of the snakeboard (radians)

𝜔𝑏

–

Angular frequency of the back platform (radian/second)

𝜔𝑓

–

Angular frequency of the front platform (radian/second)

𝜔𝑟

–

Angular frequency of the rotor (radian/second)

𝐽

–

Moment of inertia of crossbar and platforms (kg·m²)

𝑉

–

Velocity of the snakeboard‟s centre of mass (m/s)

𝑙

–

The half length of the snakeboard (m)

𝑚

–

Total mass of the system (kg)

𝑥, 𝑦 –

Coordinates of snakeboard‟s centre of mass (m)

𝛿

–

The rotor‟s angle of rotation with respect to the centre line of the snakeboard (radians)

𝜃

–

The angle between the 𝑥-axis and the central line of the snakeboard (radians)

𝜔

–

Angular frequency of platforms and rotor (radian/second)

viii

1.

Introduction 1.1.

Nonholomonic Systems

A nonholomonic system is one whose state is dependent on the trajectory (path) performed to reach that position. The parameters which control the system, described by nonholomonic constraints, may begin at a set of values, change continuously (usually over time) and return to their original values without the system having returned to its original state. A good, qualitative comparison to illustrate the difference between a nonholomonic system and a holomonic one is a sack trolley (a two wheeled device for carrying heavy loads, Figure 1) and a stationary robotic arm with two pivot joints (Figure 2).

Wheel Rotation (radians)

Right wheel Left wheel

𝑦 Initial position

𝑂

Final position

Time (seconds)

𝑥 Figure 1 – Nonholomonic system (sack trolley)

Pivot Rotation (radians)

Pivot 1 Pivot 2

𝑦 Initial position and final position

𝑂

𝑥 Figure 2 – Holomonic system (robotic arm)

–1–

Time (seconds)

In the case of the holomonic system, when the initial configuration occurs the system always returns to its original position. In contrast, when the nonholomonic system is returned to its original configuration (wheels at zero rotation) it does not have to return to its original position. Nonholomonic kinematics cannot be described by algebraic equations that relate the parameters of the system (such as wheel rotation) with the absolute position and orientation of it. Instead, differential relationships must be used. This has implications regarding system control and implementing motion planning. As the complexity of robotics increases, the desire to automate tasks normally carried out by humans grows stronger. However, from the sack trolley example, it becomes apparent that this is not a trivial task because there are an infinite number of paths the trolley can travel. Humans have an intuitive understanding of the optimal path and can react to changing circumstances. For example, if a pet walked in front of the trolley the person pushing it would automatically avoid the new obstacle and calculate a different route. This is a significant challenge for a robot and although great advances have been made, artificial intelligence in this area is below that of a human. Many robots that mimic nature (like robotic fish and snakes) are subject to nonholomonic constraints. Research has been carried out on animals with the aim of improving the motion of their robotic counterparts. Unlike many nonholomonic systems that could be turned into robotic form, the snakeboard does not have a direct biological counterpart. The author chose to analyse snakeboards because they are an interesting example of nonholomonic motion and his personal interest in the sport.

–2–

1.2.

Background to Snakeboarding

The snakeboard was invented in 1989 by James Fisher and Oliver Macleod Smith (Figure 3). The popularity of snakeboarding (the sport of riding a snakeboard) peaked in the 1990‟s, however, it has continued in different forms and names such as “streetboarding”

Figure 3 – Image from Patent Number: 4955626 relating to snakeboard

At first glance an observer would note the similarity between the snakeboard and a regular skateboard, however, the snakeboard moves with a strange snake-like motion (hence the name). Other differences also become apparent, for example, the rider need not put his foot down on the ground to propel himself. The author‟s interest was piqued at an early age by the sight of riders seemingly defying gravity by riding uphill without pushing their feet against the ground. Instead of pushing, a snakeboarder will gain momentum by coupling the twisting of the torso and rotations of the feet. This is made possible by the special construction of the snakeboard. Similar in size to the skateboard, it also has two small wheels at either end. The key difference is that the platforms upon which the rider stand are connected to the rigid crossbar with revolute joints that allow the trucks (the axle units where the wheels are attached) to rotate through the vertical axis. Figure 4 is a labelled photograph of a snakeboard and the features previously described.

–3–

Crossbar

Wheel

Revolute joint

Trucks

Platform

Figure 4 – Labelled photograph of a snakeboard

The beginner will struggle to attain any movement along the ground because the motions required of the body and feet are unintuitive. To move the torso must be rotated back and forth while simultaneously turning their feet so they point in (Figure 5) and then out (Figure 6).

Figure 5 – Feet pointing in

Figure 6 – Feet pointing out

–4–

2.

Literary Review

Academic interest in snakeboards began in 1995 with a mathematical model created by Andrew Lewis and James Ostrowski et al [1]. Various simplifications were made, including the assumption that the platforms rotate in equal and opposite directions. They also demonstrated the three basic gaits analysed in this paper. Prior to this, it was recognised that sinusoids could be used to control nonholomonic systems and Richard Murray and S. Shankar investigated this in 1993 [2]. This proved to be a starting point for others such as Stefano Iannitti and Kevin Lynch to develop motion planning methods for the snakeboard [3]. An analytical method for motion planning of a snakeboard was presented by Andrew Lewis and F. Bulllo in 2004 [4]. Other academic work regarding snakeboards includes that of A. Asnafi and M Mahzoon to develop interesting and flower-like gaits for the snakeboard [5]. Various attempts have been made to create a robotic snakeboard but the author could not find any well-documented examples except that of undergraduate electrical and control engineering student, Eddy Veltman [6]. A.S Kuleshov improved the mathematical model proposed in [1] by: “[taking] into account an opportunity that platforms of a snakeboard can rotate independently from each other” In his paper [7], he creates equations of motion derived in the Gibbs-Appell form that describe the motion of the snakeboard.

–5–

3.

Mathematical Model

To investigate the dynamics of snakeboard motion a mathematical model was created. Figure 7 represents how the snakeboard is modelled in this paper. It is assumed that it can only move in the 𝑥𝑦 plane and that 𝑂𝑥𝑦 is a fixed coordinate system. Point 𝐴, the centre of mass of the system, is given coordinates 𝑥 and 𝑦. The angle between the 𝑥-axis and the central line of the snakeboard is denoted by 𝜃. The front and back platforms can rotate independently and their angles relative to the central line of the snakeboard are denoted by 𝜑𝑓 and 𝜑𝑏 respectively. The rider is simplified to a rotor in the centre of the snakeboard. The rotor‟s angle of rotation with respect to the centre line is denoted by 𝛿.

𝜑𝑓

𝛿

𝐴

𝜃

𝜑𝑏 𝑦

𝑂

𝑥 Figure 7 – Mathematical representation

3.1.

Nonholomonic Constraints

It is assumed that the wheels may only roll without slipping on the 𝑥𝑦 plane. This restriction is a nonholomonic constraint and takes the form of Equation 1 for the front platform and Equation 2 for the back platform. 0 = − sin 𝜑𝑓 + 𝜃 𝑥 + cos 𝜑𝑓 + 𝜃 𝑦 + 𝑙 𝑐𝑜𝑠(𝜑𝑓 )𝜃

Equation 1

0 = − sin 𝜑𝑏 + 𝜃 𝑥 + cos 𝜑𝑏 + 𝜃 𝑦 − 𝑙 𝑐𝑜𝑠(𝜑𝑏 )𝜃

Equation 2

–6–

3.2.

Controlled Parameters

There are three controlled parameters in this system: 𝜑𝑓 , 𝜑𝑏 and 𝛿. By performing periodic loops with the controlled parameters, different gaits can be achieved. Gaits are defined as “a generator of motion in a certain direction”[1] . To perform these periodic loops three equations as functions of time were created for 𝜑𝑓 , 𝜑𝑏 and 𝛿 (Equation 3, Equation 4 and Equation 5 respectively): 𝜑𝑓 𝑡 = 𝑎𝑓 sin(𝜔𝑓 𝑡 + 𝛽𝑓 )

Equation 3

𝜑𝑏 𝑡 = 𝑎𝑏 sin(𝜔𝑏 𝑡 + 𝛽𝑏 )

Equation 4

𝛿 𝑡 = 𝑎𝑟 sin(𝜔𝑟 𝑡 + 𝛽𝑟 )

Equation 5

The control parameters 𝑎𝑓 , 𝑎𝑏 and 𝑎𝑟 determine the amplitude of rotation whilst the frequencies are determined by 𝜔𝑓 , 𝜔𝑏 and 𝜔𝑟 . Each controlled parameter can be shifted in phase using the respective control parameter 𝛽𝑓 , 𝛽𝑏 and 𝛽𝑟 . The periodic loops are used in all sections of this paper except the physical input controller where the controlled parameters are taken directly.

3.3.

Physical Parameters

Table 1 lists the values of parameters used in this investigation. They represent the physical properties of a snakeboard. These were constant throughout the paper for every simulation. Parameter

Value

l m J Jr Jp

0.285 m 75 kg 0.22 kg m2 14 kg m2 0.013 kg m2

Table 1 – Physical Parameters

3.4.

ODE45 Differential Solver

The equations of motion for the snakeboard are ordinary differential equations. To solve these, the MATLAB ODE45 function was used. The ODE45 solver is recommended for nonstiff problem types and has a good –7–

order of accuracy. It requires the initial conditions, for this problem 𝑥, 𝑦, 𝜃 and 𝑉, in vector form. There are two outputs, a column vector of time points and a solution array with rows corresponding to the solution at the relevant time point. For this investigation, the position and velocity of the snakeboard between the initial and final solution are of interest. This is because the motion is unintuitive so it is beneficial for the interim positions to be known. Because the ODE45 solver is an adaptive integrator, only an initial time and final time are required. However, to obtain an output that describes the full motion a timestep can be applied. Small timesteps will give smooth trajectories in contrast to large timesteps where the trajectory of the snakeboard is unclear. It should be emphasised that decreasing the timestep has no appreciable effect on the accuracy of the final solution (see Table 2) and the computational time is only increased because more timesteps need to be stored in the computer‟s memory. A test was conducted were everything remained constant except timestep with the results shown in Table 2. The final 𝑥 position was recorded to prove accuracy was unaffected and the solution time for the integration was measured to check if there were any significant differences in computational time. Timestep (s) ODE Solution time (s) Final 𝒙 position (m)

5 2.5 1 0.5 0.25 0.1 0.05 0.01

0.3183 0.3373 0.3329 0.3339 0.3357 0.3470 0.3569 0.4085

0.9608 0.9608 0.9608 0.9609 0.9608 0.9609 0.9608 0.9608

Table 2 – ODE 45 Timestep

Table 2 confirms that there is no increase in accuracy when time step is reduced but the solution time increases. Therefore, a qualitative approach was taken. Trajectory plots for each timestep were created (see Appendix I – Timestep) and it was decided that 0.1s was the largest possible timestep to give satisfactory results. –8–

3.5.

Dynamical Equations of Motion

The system of equations derived for the motion of a snakeboard derived in [7] was rearranged into four, coupled, non-linear differential equations (Equation 6, Equation 7, Equation 8 and Equation 9) that were in a form suitable for MATLAB to solve. 𝑥 = 𝑉 cos 𝜃 −

𝑉 sin 𝜓2 𝑠𝑖𝑛 𝜃 𝑐𝑜𝑠𝜓1 + 𝑐𝑜𝑠𝜓2

Equation 6

𝑦 = 𝑉 sin 𝜃 +

𝑉 sin 𝜓2 𝑐𝑜𝑠 𝜃 𝑐𝑜𝑠𝜓1 + 𝑐𝑜𝑠𝜓2

Equation 7

𝜃=

− 𝑉=

𝑉 sin 𝜓1 𝑙 (cos 𝜓1 + cos 𝜓2 )

(𝑑1 𝛿 + 𝑑2 𝜓2 ) sin 𝜓1 − 𝑉 × 𝑀(𝑡) (cos 𝜓1 + cos 𝜓2 ) sin2 𝜓2 + 𝑘 2 sin2 𝜓1 1+ cos 𝜓1 + cos 𝜓2 2

Equation 8

Equation 9

Where: 𝑀 𝑡 =

𝜓2 sin 𝜓2 cos 𝜓2 + 𝑘 2 𝜓1 sin 𝜓1 cos 𝜓1 cos 𝜓1 + cos 𝜓2 2 +

𝜓1 sin 𝜓1 𝜓2 sin 𝜓2 (sin2 𝜓2 + 𝑘 2 sin2 𝜓1 ) cos 𝜓1 + cos 𝜓2 3

Equation 10

𝐽 + 𝐽𝑟 + 2𝐽𝑝 𝑚𝑙 2

Equation 11

𝑑1 =

𝐽𝑟 𝑚𝑙

Equation 12

𝑑2 =

𝐽𝑝 𝑚𝑙

Equation 13

𝑘2 =

𝜓1 = 𝜑𝑓 − 𝜑𝑏

Equation 14

𝜓2 = 𝜑𝑓 + 𝜑𝑏

Equation 15

–9–

3.6.

Validation

To validate the mathematical model had been correctly entered into MATLAB two of the results from [7] were replicated. The method of applying controlled parameters is slightly different from the one used in this investigation. Instead of having the individual functions for each platform with respect to time, they are merged together (Equation 14 and Equation 15). This method is limited because it is harder to define the periodic controls. Most of the other papers use a similar method to the one used in the rest of this paper. 3.6.1.

Validation Test 1

This case was labelled as an example of a “nonresonant case” which is representative of beginners unable to ride a snakeboard. There is only a small displacement in the 𝑥-direction and no appreciable forward motion. The controlled parameters are given in Table 3 and the trajectory calculated by the MATLAB model used for this paper is shown in Figure 9. Figure 8 is the trajectory plotted in [7]. The trajectories are similar, both quantitatively and qualitatively. Controlled Parameter

Function

0.3 sin

1

𝑡 3 1 0.2 sin 𝑡 2 1 0.7 sin 𝑡 2

𝜓1 𝑡 𝜓2 (𝑡) 𝛿(𝑡)

Table 3 – Validation Test 1 parameters

Figure 8 – Trajectory of snakeboard in Validation Test 1 from [7]

– 10 –

Figure 9 – Trajectory of snakeboard in Validation Test 1 using MATLAB model in this paper

3.6.2.

Validation Test 2

This test was similar to the previous one but used a resonant gait. The controlled parameters are given in Table 4 and the trajectory calculated by the MATLAB model used for this paper is shown in Figure 10. Figure 11 is the trajectory plotted in [7]. Controlled Parameter

Function

0.3 sin

1

𝑡 3 1 0.2 sin 𝑡 2 1 0.7 sin 𝑡 2

𝜓1 𝑡 𝜓2 (𝑡) 𝛿(𝑡)

Table 4 – Validation Test 1 parameters

Figure 10 – Trajectory of snakeboard in Validation Test 1 using MATLAB model in this paper

– 11 –

Figure 11 – Trajectory of snakeboard in Validation Test 1 from [7]

Again, the trajectories are both similar quantitatively and qualitatively. From these two tests, it was concluded that the MATLAB model used in this paper replicated the results of previous work sufficiently well and validated the model.

3.7.

Limitations

The major limitation of the dynamic equations used in this paper is friction is not taken into account. From the author‟s experience of riding a snakeboard, this is unrealistic, especially at lower speeds. An example of the lack of realism from neglecting friction is when [7] concludes that: Even if the rider doesn’t rotate his torso, he can propel [the] snakeboard forward using only his legs This may be true in an idealised world but an engineer designing for practical applications would be unable to achieve motion without a rotor of a reasonable moment of inertia. An experienced rider can shift their body weight allowing for tight turns at high speeds, however, there is a limit to what is possible without the rider falling off. This simulation assumes that the rider cannot fall off or the robot remains stable at all times. If a robotic snakeboard was being designed then the centre of mass of the rotor should be as low as possible to minimise the chance of tipping. This model also assumes that the surface the snakeboard is perfectly smooth and level. In reality, a gradient always exists and the surface undulates. Future work could take into account these effects but it would add a significant level of complexity onto the model. One phenomenon observed by the author when riding downhill at high speeds is “speed wobble”. The snakeboard platforms begin to oscillate rapidly and uncontrollably. However, this is outside the scope of this model and study.

– 12 –

3.8.

Plotting and Visualisation

The majority of the figures in this paper are line graphs because they are the most suitable. However, other methods were created to represent the motion of the snakeboard over time. These were useful because the snakeboard‟s motion can be unintuitive so it is often difficult to understand a line plot of its trajectory. Static Overlay

3.8.1.

An example of the plot defined as “Static Overlay” for the purposes of this paper is shown in Figure 12. Two-dimensional snapshots of the snakeboard‟s position and state are taken at fixed intervals and overlaid to form a single image.

Figure 12 – Static Overlay visualisation example

An indication of the snakeboard‟s velocity can be observed from the separation between each snapshot; a large displacement corresponds to a large velocity. The main problem with this visualisation is when snapshots overlap, so distinguishing between them becomes difficult.

– 13 –

3.8.2.

Animated 2D

This visualisation is similar to the previous one but instead of overlaying each snapshot, they are shown sequentially. This makes it unsuitable for a printed medium but an example is shown in Figure 13.

Figure 13 – Animated 2D visualisation example

3.8.3.

Animated 3D

This visualisation is a 3D version of the “Animated 2D” version. An example is shown in Figure 14. The camera position could be fixed in space or move relative to the snakeboard‟s position and angle.

Figure 14 – Animated 3D visualisation example

3.9.

MATLAB Code

The MATLAB code written by the author to produce the results in this paper can be found in “Appendix II – MATLAB Code.”

– 14 –

4.

Resonance Gait

Resonance is defined as when the platforms and rotor oscillate at the same frequency, that is, when 𝜔𝑓 , 𝜔𝑏 and 𝜔𝑟 are at a ratio of 1:1:1 and there is no phase shift (𝛽𝑓 = 𝛽𝑏 = 𝛽𝑟 = 0). It is the gait used by snakeboarders and produces a net displacement in the 𝑥-direction. The amplitudes 𝑎𝑓 and 𝑎𝑏 must be in opposite directions (but not necessarily equal magnitude) for motion to occur. As an example of resonance gait, Figure 15 shows the trajectory of the snakeboard when the simulation is run with the control parameters listed in Table 5. The controlled parameters over time are shown in Figure 16. Parameter

Value

𝑎𝑓 𝑎𝑏 𝑎𝑟 𝜔𝑓 𝜔𝑏 𝜔𝑟 𝛽𝑓 𝛽𝑏 𝛽𝑟

0.3 rad -0.3 rad 0.7 rad 1 rad/s 1 rad/s 1 rad/s 0 rad 0 rad 0 rad

Table 5 – An example of resonance gait control parameters.

Figure 15 – Example of resonance trajectory

Figure 16 – Controlled parameters during resonance

– 15 –

Direction of Travel

4.1.

The parameters in Table 6 were kept constant while various combinations of platform amplitudes were tested. The results of the tests are shown in Table 7. It can been seen that when the direction of 𝑎𝑓 and 𝑎𝑟 are the same the snakeboard accelerate in the direction of the front platform. Likewise, when 𝑎𝑏 and 𝑎𝑟 are the same the snakeboard will accelerate in the direction of the back platform. Parameter

Value

𝜔𝑓 𝜔𝑏 𝜔𝑟 𝛽𝑓 𝛽𝑏 𝛽𝑟

1 rad/s 1 rad/s 1 rad/s 0 rad 0 rad 0 rad

Table 6 – Direction of travel during resonance gait, constant parameters

𝒂𝒇 (rad)

𝒂𝒃 (rad)

0.3 0.3 -0.3 -0.3

-0.3 -0.3 0.3 0.3

𝒂𝒓 (rad) 𝒙-displacement Trajectory

-0.7 0.7 0.7 -0.7

Negative Positive Negative Positive

Figure 17 Figure 18 Figure 19 Figure 20

Table 7 – Direction of travel for resonance gait

Figure 17

Figure 18

Figure 19

Figure 20

– 16 –

4.2.

Amplitude

To gain an understanding of how the amplitude of the rotor and platforms affected the motion of the snakeboard, surface plots were created. The distance travelled in the 𝑥-direction was calculated for various time periods. For each calculation, the platform amplitudes were assumed to have equal magnitudes in opposite directions (𝑎𝑓 = −𝑎𝑏 ). The platform and rotor amplitudes were varied from 0 to 1.5 radians in increments of 0.01 radians. The maximum limit of 1.5 radians was chosen because rotation greater than 90 degrees in the platforms will cause the snakeboard to exhibits strange behaviour and the gait is not resonant. The parameters in Table 8 were kept constant and the results are listed in Table 9. Parameter

Value

𝜔𝑓 𝜔𝑏 𝜔𝑟 𝛽𝑓 𝛽𝑏 𝛽𝑟

1 rad/s 1 rad/s 1 rad/s 0 rad 0 rad 0 rad

Table 8 – Constant parameters during resonance amplitude tests

Simulation Time (s)

Platform amplitude associated with largest 𝒙displacement

Rotor amplitude associated with largest 𝒙displacement

Largest 𝒙displacement

Surface Plot

10 30 70 120

0.53 0.50 0.37 0.29

1.5 1.5 1.5 1.5

0.69 5.80 25.4 60.27

Figure 21 Figure 22 Figure 23 Figure 24

Table 9 – Amplitude surface plot results

Figure 21 – Amplitude surface plot for 10 seconds

Figure 22 – Amplitude surface plot for 30 seconds

– 17 –

Figure 23 – Amplitude surface plot for 70 seconds

Figure 24 – Amplitude surface plot for 120 seconds

From Table 9 it can be seen that the greater the largest rotor amplitude always gives the largest 𝑥-displacement. However, when starting from rest, if large platform amplitudes are used for long times then the oscillations of the snakeboard in the 𝑦-direction grow very large. This is best explained by the trajectory shown in Figure 25. The large oscillations in the 𝑦-direction soon become unstable and the snakeboard doubles back on itself meaning the displacement in the 𝑥-direction is reduced. The concept of stability in the resonance gait is investigated in section 4.5 titled “Instability”

Figure 25 – Trajectory during resonance gait with large platform amplitudes

For short run times a large platform amplitude is desirable because it gets the snakeboard moving. Optimising the platform amplitudes to maximise displacement in the 𝑥-direction is investigated further in section 11 titled “Stable Resonance Gait”.

– 18 –

4.3.

Angular Frequency

The effect of changing the angular frequency of the platforms and rotor was investigated in a similar way to that of the amplitude. The angular frequencies 𝜔𝑓 , 𝜔𝑏 and 𝜔𝑟 were varied from 0 to 2 rad/s in increments of 0.05 rad/s. The parameters in Table 10 were kept constant and a run time of 200 seconds for each combination was used. A surface plot of the results is shown in Figure 26. Parameter

Value

𝑎𝑓 𝑎𝑏 𝑎𝑟 𝛽𝑓 𝛽𝑏 𝛽𝑟

0.3 rad -0.3 rad 0.7 rad 0 rad 0 rad 0 rad

Table 10 – Constant parameters during resonance gait, angular frequency test

Figure 26 – Surface plot of angular velocities

This plot shows that there is no other combination of angular frequencies apart from resonance that produces a significant net displacement in the 𝑥direction. It also indicates how sensitive resonance is. If the ratio of platforms and rotors is not extremely close to 1:1:1 then its associated gait is not achieved. The small “humps” that are visible occur when 𝜔𝑓 , 𝜔𝑏 and 𝜔𝑟 are at a ratio of 1:1:2. This ratio produces the rotation gait and is discussed in section 5.

– 19 –

Assuming the ratio of 𝜔𝑓 , 𝜔𝑏 and 𝜔𝑟 is 1:1:1 then the greater they are the further the snakeboard will travel in a given period until instability occurs (Figure 27).

Figure 27 – Displacement during a range of angular frequencies

There is an exponential relationship between angular frequency and the displacement travelled in the 𝑥-direction until 𝜔 is greater than 2.25 rad/s in this case. The reason it does not continue to grow is that the gait becomes unstable. For motion planning care needs to be taken that if the angular frequencies were increased to travel that the gait does not become unstable if long runtimes are used.

4.4.

Phase shift

Another surface plot (Figure 28) was created to investigate the effect of changing the phase shift of the platforms and rotors. The phase shift of the rotor and the back platform (𝛽𝑟 and 𝛽𝑏 ) was varied from 0 to 4𝜋 radians in increments of 0.05 radians. The simulation run time was 50 seconds for each combination and the parameters in Table 11 were kept constant. Parameter

Value

𝑎𝑓 𝑎𝑏 𝑎𝑟 𝜔𝑓 𝜔𝑏 𝜔𝑟

0.3 rad -0.3 rad 0.7 rad 1 rad/s 1 rad/s 1 rad/s

Table 11 – Constant parameters for phase shift test

– 20 –

Figure 28 – Surface plot of the effects of phase shift

The surface plot confirms the conclusions made about the direction of motion. A phase shift of 2𝜋 radians is equivalent to reversing the direction of the amplitudes.

4.5.

Instability

For simulations with long run times, the snakeboard‟s trajectory becomes unstable. This was noted in [1]: “one cannot just increase the run time and input magnitude to get longer distances traveled in the (1,1,1) [resonance] gait” Figure 29 shows the trajectory of the snakeboard during a resonance simulation with a long run time and Figure 30 is a close-up of the unstable part. Unlike the previous example of instability (Figure 25), more reasonable platform amplitudes were used (𝑎𝑓 = −𝑎𝑏 = 0.3) and instability still occurs, albeit longer into the run time. Although efforts were made to determine a method of predicting when instability would occur this was not achieved. Instead, checks can be made after the trajectory has been calculated to ensure the snakeboard has not travelled backwards.

– 21 –

Figure 29 – Trajectory becoming unstable

Figure 30 – Close-up of where trajectory becomes unstable

– 22 –

5.

Rotation Gait

The rotation gait was discovered first in [1] and occurs when the angular frequencies 𝜔𝑓 , 𝜔𝑏 and 𝜔𝑟 are at a ratio of 1:1:2. It generates a net motion in the 𝜃-direction and when combined with an initial non-zero velocity, circular motion is produced. Snakeboard riders do not use this method to turn; instead, they slightly adapt the resonance gait to correct the direction they wish to go. Nonetheless, robotic snakeboards would not struggle to use this gait so it was analysed. As an example of the rotation gait, Figure 45 shows the trajectory of the snakeboard when the simulation is run with the control parameters listed in Table 12. The controlled parameters over time are shown in Figure 45.

Parameter

Value

𝑎𝑓 𝑎𝑏 𝑎𝑟 𝜔𝑓 𝜔𝑏 𝜔𝑟 𝛽𝑓 𝛽𝑏 𝛽𝑟

0.3 rad -0.3 rad 0.7 rad 1 rad/s 1 rad/s 2 rad/s 0 rad 0 rad 0 rad

Table 12 – An example of resonance gait control parameters

Figure 31 – Trajectory of snakeboard in rotation gait

Figure 32 – Controlled parameters during rotation gait

– 23 –

Direction

5.1.

The first concern was to understand how to control the direction of the rotation. From the work done in section 4 on resonance it seemed likely that the directions of 𝑎𝑓 , 𝑎𝑏 and 𝑎𝑟 would influence the direction of rotation. Just as no movement occurs when the platform amplitudes are in the same direction during the resonance gait, the same is true for the rotation gait. Consequently, there were four combinations of amplitude directions and these were each tested in four separate simulations. The parameters in Table 13 were kept constant for each of the four simulations. The simulation run time was fixed at 25 seconds and the results are shown in Table 14. Parameter

Value

𝜔𝑓 𝜔𝑏 𝜔𝑟 𝛽𝑓 𝛽𝑏 𝛽𝑟

1 rad/s 1 rad/s 2 rad/s 0 rad 0 rad 0 rad

Table 13 – Constant parameters during rotation direction test

𝒂𝒇 (rad)

𝒂𝒃 (rad)

𝒂𝒓 (rad)

Rotation Direction

Trajectory

0.3 0.3 -0.3 -0.3

-0.3 -0.3 0.3 0.3

-0.7 0.7 0.7 -0.7

Clockwise Counter-clockwise Counter-clockwise Clockwise

Figure 33 Figure 34 Figure 35 Figure 36

Table 14 – Results of rotation direction test

Figure 33

Figure 34

– 24 –

Figure 35

Figure 36

The results in Table 14 show that as long as the amplitudes of the platforms are opposite in direction, it is the direction 𝑎𝑟 that determines the direction of rotation. A positive value for 𝑎𝑟 causes a counter-clockwise rotation and a negative value for 𝑎𝑟 causes clockwise rotation.

Rate of rotation

5.2.

The rate of rotation can be changed by varying the amplitudes or angular frequencies of the controlled parameters. Both of these were investigated in the following sections. 5.2.1.

Angular frequency

To investigate the effect of angular frequency a script was created that varied 𝜔𝑓 and 𝜔𝑏 from 0 to 2 in increments of 0.1 radians while keeping a ratio of: 𝜔𝑟 = 2 𝜔𝑓 = 2 𝜔𝑓 . Each iteration was run for 25 seconds and the angle of rotation, 𝜃 was recorded and plotted (Figure 37). The parameters in Table 15 were kept constant. It can be seen that the relationship between rate of angular rotation and the angular frequencies is linear.

Parameter

Value

𝑎𝑓 𝑎𝑏 𝑎𝑟 𝛽𝑓 𝛽𝑏 𝛽𝑟

0.3 rad -0.3 rad -0.3 rad 0 rad 0 rad 0 rad

Table 15 – Control parameters for rate of rotation, angular frequency test

Figure 37 – Rate of rotation for a range of angular frequencies

– 25 –

5.2.2.

Amplitude of platforms and rotors

The first amplitude investigation was changing 𝑎𝑟 from 0 to 0.7 radians in increments of 0.01 radians. Each iteration was run for 25 seconds and the angle of rotation, 𝜃 was recorded and plotted (Figure 38). The parameters in Table 16 were kept constant. The relationship between the rate of angular rotation and amplitude of the rotor rotation is linear. Parameter

Value

𝑎𝑓 𝑎𝑏 𝜔𝑓 𝜔𝑏 𝜔𝑟 𝛽𝑓 𝛽𝑏 𝛽𝑟

0.3 rad -0.3 rad 1 rad/s 1 rad/s 2 rad/s 0 rad 0 rad 0 rad

Table 16 – Constant parameters for rate of rotation rotor amplitude test

Figure 38 – Rate of rotation for a range of rotor amplitudes

The second amplitude investigation was changing 𝑎𝑓 and 𝑎𝑏 from 0 to 0.7 radians in increments of 0.01 radians. Each iteration was run for 25 seconds and the angle of rotation, 𝜃 was recorded and plotted (Figure 39). The parameters in Table 17 were kept constant. It can be seen that the relationship between the rate of angular rotation and amplitude of the rotor rotation is non-linear.

Parameter

Value

𝑎𝑟 𝜔𝑓 𝜔𝑏 𝜔𝑟 𝛽𝑓 𝛽𝑏 𝛽𝑟

0.7 rad 1 rad/s 1 rad/s 2 rad/s 0 rad 0 rad 0 rad

Table 17 – Constant parameters for rate of rotation platform amplitude test

Figure 39 – Rate of rotation for a range of platform amplitudes

– 26 –

5.2.3.

Initial Velocity

Rather than starting with an initial velocity of zero, the parking gait was investigated with an initial velocity. The parameters in Table 18 were used and an initial velocity of 0.5 m/s was applied. The simulation was run for 60 seconds and the trajectory of the snakeboard is shown in Figure 40.

Parameter

Value

𝑎𝑓 𝑎𝑏 𝑎𝑟 𝜔𝑓 𝜔𝑏 𝜔𝑟 𝛽𝑓 𝛽𝑏 𝛽𝑟

0.3 rad -0.3 rad 0.7 rad 1 rad/s 1 rad/s 2 rad/s 0 rad 0 rad 0 rad

Table 18 – Constant parameters for initial velocity test

Figure 40 – Trajectory with an initial velocity applied

Then an identical simulation was run except the angular frequencies were increased so the parameters in Table 19 were used. The trajectory is shown in Figure 41.

Parameter

Value

𝑎𝑓 𝑎𝑏 𝑎𝑟 𝜔𝑓 𝜔𝑏 𝜔𝑟 𝛽𝑓 𝛽𝑏 𝛽𝑟

0.3 rad -0.3 rad 0.7 rad 1.5 rad/s 1.5 rad/s 3 rad/s 0 rad 0 rad 0 rad

Table 19 – Constant parameters for initial velocity test

Figure 41 – Trajectory with an initial velocity applied

It can be seen that increasing the angular frequency reduces the radius of the circular trajectory. Attempts were made to quantify this change but without success. This is an area for further work because an understanding of it and the resonant gait would allow full motion planning. – 27 –

6.

Parking Gait

The “parking” gait is the third and final one discussed in [1]. It occurs when the angular frequencies 𝜔𝑓 , 𝜔𝑏 and 𝜔𝑟 are at a ratio of 2:2:3. It produces a net displacement in the 𝑦-direction. Like the rotation gait, the parking gait is not one used by snakeboard riders. As an example of the parking gait, Figure 42 shows the trajectory of the snakeboard when the simulation is run with the control parameters listed in Figure 21. The controlled parameters over time are shown in Figure 43.

Parameter

Value

𝑎𝑓 𝑎𝑏 𝑎𝑟 𝜔𝑓 𝜔𝑏 𝜔𝑟 𝛽𝑓 𝛽𝑏 𝛽𝑟

0.3 rad -0.3 rad 0.7 rad 2 rad/s 2 rad/s 3 rad/s 0 rad 0 rad 0 rad

Table 20 – Constant parameters for parking gait example

Figure 42 – Trajectory of snakeboard during parking gait

Figure 43 – Control Parameters during parking gate

6.1.

Direction

Like the other gaits the first priority was to understand how to control the direction of net displacement in the 𝑦-axis. The four combinations of amplitude values which produced at net displacement were tested in a

– 28 –

simulation with a run time of 15 seconds (Table 22). The parameters in Table 21 were kept constant. Parameter

Value

𝜔𝑓 𝜔𝑏 𝜔𝑟 𝛽𝑓 𝛽𝑏 𝛽𝑟

2 rad/s 2 rad/s 3 rad/s 0 rad 0 rad 0 rad

Table 21 – Constant parameters for parking gait direction test

𝒂𝒇 (rad)

𝒂𝒃 (rad)

0.3 0.3 -0.3 -0.3

-0.3 -0.3 0.3 0.3

𝒂𝒓 (rad) Net displacement in the 𝒚-direction Trajectory

-0.7 0.7 0.7 -0.7

Positive Positive Negative Negative

Figure 44 Figure 45 Figure 46 Figure 47

Table 22 – Results for parking gait direction test

Figure 44

Figure 45

Figure 46

Figure 47

– 29 –

Unlike the rotation and resonance gait there is no simple rule for determining which direction the net 𝑦-displacement will be in. Also, each of the possible directions for the parking gait can be offset in the positive or negative 𝑥direction.

6.2.

Rate of change in the 𝒚-direction

The rate of change in the 𝑦-direction can be changed by varying the amplitudes or angular frequencies of the controlled parameters, and are investigated in the following sections. 6.2.1.

Angular Frequency

To investigate the effect of angular frequency a script was created that varied 𝜔𝑓 and 𝜔𝑏 from 0 to 2 in increments of 0.1 radians while keeping a ratio of 3𝜔𝑟 = 2 𝜔𝑓 = 2 𝜔𝑓 . Each iteration was run for 25 seconds and the 𝑦displacement was recorded and plotted (Figure 48). The parameters in Table 23 were kept constant. The relationship between rate of change in the 𝑦direction and the angular frequencies is reasonably linear but unstable in the larger angular frequencies.

Parameter

Value

𝑎𝑓 𝑎𝑏 𝑎𝑟 𝛽𝑓 𝛽𝑏 𝛽𝑟

0.3 rad -0.3 rad -0.7 rad 0 rad 0 rad 0 rad

Table 23 – Constant parameters for angular frequency variation

Figure 48 – Displacement for a range of platform frequencies

– 30 –

6.2.2.

Amplitude of platforms and rotors

The first amplitude investigation was changing 𝑎𝑟 from 0 to 0.7 radians in increments of 0.01 radians. Each iteration was run for 25 seconds and the 𝑦displacement was recorded and plotted (Figure 49). The parameters in Table 24 were kept constant. The investigation shows that the relationship between the rate of angular rotation and amplitude of the rotor rotation is non-linear. Parameter

Value

𝑎𝑓 𝑎𝑏 𝜔𝑓 𝜔𝑏 𝜔𝑟 𝛽𝑓 𝛽𝑏 𝛽𝑟

0.3 rad -0.3 rad 1 rad/s 1 rad/s 2 rad/s 0 rad 0 rad 0 rad

Table 24 – Constant parameters for rotor amplitude variation

Figure 49 – Displacement for a range of rotor amplitudes

The second amplitude investigation was changing 𝑎𝑓 and 𝑎𝑏 from 0 to 0.7 radians in increments of 0.01 radians. Each iteration was run for 25 seconds and the 𝑦-displacement was recorded and plotted (Figure 50). The parameters in Table 25 were kept constant. The investigation shows that the relationship between the rate of angular rotation and amplitude of the rotor rotation is non-linear.

Parameter

Value

𝑎𝑟 𝜔𝑓 𝜔𝑏 𝜔𝑟 𝛽𝑓 𝛽𝑏 𝛽𝑟

0.7 rad 1 rad/s 1 rad/s 2 rad/s 0 rad 0 rad 0 rad

Table 25 – Constant parameters for platform amplitude variation

Figure 50 – Displacement for a range of platform amplitudes

– 31 –

6.3.

Initial Velocity

Rather than starting with an initial velocity of zero, the parking gait was investigated with an initial velocity. The parameters in Table 20 were used and an initial velocity of 0.5 m/s was applied. The simulation was run for 100 seconds and the trajectory of the snakeboard is shown in Figure 51. The velocity over time is plotted in Figure 52.

Figure 51 – Trajectory of snakeboard when an initial velocity is applied to the parking gait

Figure 52 – Snakeboard velocity when an initial velocity is applied to the parking gate

Once an initial velocity is applied, the trajectory is in a straight line, similar to resonance gait. Interestingly the average velocity stays relatively constant over time unlike the resonance gait.

– 32 –

7.

Direct Simulation Controller

Most of the research on snakeboard dynamics has focused on using mathematical tools to analyse how they behave. Using sinusoids to represent the periodic motions came from observation but to date no serious attempt at recording a real snakeboard has been made. Ideally, the inputs, platform angle and body position would be recorded as a function of time while simultaneously recording the trajectory of the snakeboard. From this, an understanding of how a real rider negotiates obstacles, generates momentum and controls the snakeboard could be gained. Applying this knowledge back into computer simulations could lead to better motion planning algorithms and further the field of nonholomonic robots, like the bio-inspired ones mentioned in the introduction. For this investigation, various methods of recording the movements of a real snakeboard were proposed and evaluated. The most accurate method would be using motion capture technology like that used in the film industry. There are facilities at the University of Strathclyde to do this, however, it was decided that it was not suitable for this project. Getting the data recorded into a form that could be analysed would be a time consuming and difficult task. In addition, access to the equipment would be limited. Another option was to attach recording devices to a snakeboard. The problem with this approach is finding a way to record the path. Accurately determining the absolute position of the snakeboard is not a trivial task. The method chosen was to construct a small device with three rotational inputs (front platform, back platform and body/rotor rotation). This device was connected to a computer running with a MATLAB simulation that took the inputs and in real-time used them in the mathematical model of the snakeboard. A graphical output was shown to user (Figure 14). It was found that having the camera move relative to the snakeboard in the 𝑥 and 𝑦 direction was beneficial. It was possible to control the simulation when the camera rotated relative to 𝜃, however, it was distracting so this was not taken forward.

– 33 –

Another benefit of the construction this device is the possibility for future work. It is ideally suited to directly controlling a robotic snakeboard and testing its behaviour. Standard input devices such as joysticks, keyboards and mice would not be suitable and limited to a method similar to that used in “Open loop periodic controls”. Get potentiometer values and convert to radians for each controlled parameter

Get the time period since timer was started . Reset the timer and start it

Update the visual display of the snakeboard

Solve the equations of motion to obtain the path of the snakeboard for the time period

Use the time period, previous controlled parameters and the new ones to calculate velocity and acceleration of controlled parameters

Figure 53 – Direct Simulation Controller Flowchart

7.1.

Limitations

It is assumed that the mathematical model is suitably accurate. If there are errors with the model, the user might unconsciously compensate for them. However, any obvious errors would be apparent to someone familiar to riding a real snakeboard so this controller acts as a validation to the mathematical model. An obvious difference between this input device and a real snakeboard is that hands are used instead of standing on a board. It was assumed that this should not be an insurmountable problem and in practice this proved to be the case. After some practicing, it became possible for the author to use his hands to control the virtual snakeboard with reasonable proficiency. It was observed that other users struggled to produce any kind of motion, just like – 34 –

beginner snakeboard riders. This is probably because the experience from riding a real snakeboard translates to controlling a virtual one. The lack of feedback apart from visual (the animation on screen) is a limitation because a rider‟s inputs are influenced by the forces felt by their body moving.

7.2.

Construction

7.2.1.

Components

Table 26 lists the components used to create the controller and their cost. The overall total was £8.28. Parts were chosen for their low cost and availability. Component

Quantity Unit Cost (£) Cost (£)

PICAXE 08M Linear Potentiometer Potentiometer Shaft Stripboard 10k Resistor 22k Resistor 3.5mm Socket LED Enclosure Mass Wire Wooden Lollipops Total Cost

1 3 3 1 1 1 1 1 1 1 4

1.50 0.50 0.80 1.30 0.10 0.10 0.08 0.10 2.00 -

1.50 1.50 1.60 1.30 0.10 0.10 0.08 0.10 2.00 £8.28

Table 26 - Components

7.2.2.

Microcontroller

The purpose of the microcontroller is to convert the analogue signals from the potentiometers into a digital signal that is transmitted to MATLAB through a serial connection. The PICAXE 08M microcontroller was chosen because of its low cost and is easily programmed. 7.2.3.

Final Design

A labelled photograph of the final design is shown in Figure 54.

– 35 –

Serial connection to Batteries computer

Mass

Platform potentiometer Rotor potentiometer Platform potentiometer LED

Circuit board

Figure 54 – Labelled photograph of the physical input controller

7.2.4.

Improvements

The potentiometers are not designed for continual use as is required for this device. Over time, their accuracy will decrease until the controller no longer functions correctly. It was decided that this was not an issue for this project because of its limited duration. However, if intensive long term research were to be carried out then a different method of measuring the angular rotation would needed such as an LVDT (linear variable differential transformer) which does not use sliding contacts.

7.3.

Results

Some practice was needed before the author was able to use the controller. From observing others unfamiliar with the concept of snakeboarding, it is apparent that being able to ride a real snakeboard is of great benefit when trying to use this. Interestingly the author showed the same dominance in his hands as he did with his feet. In skateboarding terminology, a person is either “regular” or “goofy” (there is a minority equally adept to both styles just as some ambidextrous). If your left foot is leading in the direction of travel then you‟re “regular”. Likewise, if your right foot leads then you‟re “goofy”. Having learnt

– 36 –

how to control a snakeboard with your feet, the skills are probably “transferred” to your hands. Accordingly, the simulation was set up to ensure the author caused the snakeboard to move in the positive 𝑥 direction to enable a simple analysis. 7.3.1.

Validation Test

The first test was to perform a simple validation that the controller and MATLAB code was working as intended. The front platform was given a rotational displacement in an anticlockwise direction and the back platform was given a rotational displacement of similar rotational displacement in a clockwise direction. The rotor was rotated rapidly in a clockwise direction. This should cause the snakeboard to accelerate anticlockwise (conservation of angular momentum). Then the rotor was rotate slowly back to its starting position which should return the snakeboard to zero velocity. The controlled parameters as a function of time are shown in Figure 55 and the velocity is shown in Figure 56.

Figure 55 – Controlled parameters for Validation Test

– 37 –

Figure 56 – Snakeboard velocity during Validation Test

The snakeboard behaved as expected during this test, which is a good indication the simulation is working correctly. 7.3.2.

Resonance imitation test

The aim of this test is to try to imitate the resonance gait. The results shown in (Figure 57, Figure 58 and Figure 59) are representative of the typical attempts made.

Figure 57 – Controlled parameters for Resonance imitation test

– 38 –

Figure 58 – Trajectory of the snakeboard for Resonance imitation test

Figure 59 – Velocity of the snakeboard during Resonance imitation test

The trajectory of the snakeboard is clearly not a replication of the ones produced by the controlled parameters created using periodic controls. The oscillations about the 𝑦-axis vary greatly instead of gradually increasing. This can be expected as the controlled parameters created with a sinusoidal function can be precisely controlled whereas the controller requires direct human input. The sensitivity of the snakeboard to the controlled parameters is explored further in section 9 regarding bifurcation. From Figure 59 it can be seen that the snakeboard initially accelerates rapidly before settling to oscillations about 0.8m/s. Attempts to increase the velocity proved impossible. This suggests that it is relatively easy to maintain – 39 –

a certain velocity but trying to go faster is significantly more difficult. Trying to quantify this is difficult because by their nature, planned controlled parameters are consistent. Future work could try to replicate the inaccuracy of controlled parameters by applying a degree of randomness. This would be beneficial because nonholomonic systems in the real world will always have some imprecision. This could be in the physical components such as motors or control systems. Using the example of a robotic snakeboard, if the platforms could not be positioned accurately then its velocity would be limited. 7.3.3.

Rotation Imitation Test

The aim of this test is to try to imitate the rotation gait. The results shown in (Figure 60, Figure 61 and Figure 62) are representative of the typical attempts made.

Figure 60 – Controlled parameters for Rotation Imitation Test

Figure 61 – Trajectory for Rotation Imitation Test

– 40 –

Figure 62 – Velocity of the snakeboard during Rotation Imitation Test

Although rotation is achieved, the trajectory is erratic. Similarly, the controlled variables are unclear, appearing similar to those used for the resonance gait except for some variations. This makes sense, as a resonance gait combined with small alterations is how a human riding a real snakeboard would control it. The snakeboard maintains a reasonably constant average speed.

– 41 –

8.

Open loop periodic controls

Whilst experimenting with varying the controlled parameters it was noticed that trying different combinations to observe the behaviour of the snakeboard was taking a long time. Each simulation has to have its parameters adjusted in a text file before being run and the results viewed. This does allow accuracy and control which is why the majority of the work done in this paper used this method. The open loop, periodic controls developed and covered in this section sacrifices precision for ease of use. The periodic inputs for the rotor and the platforms can be changed using a GUI (graphic user interface) which controls all the parameters, shown in Figure 63.

Figure 63 – GUI for “Open loop periodic controls”

While the user changes parameters in the GUI a real-time 3D model is displayed (Figure 14). This allows the user to observe the behaviour as changes are made which facilitates experimentation. This method is much easier to use than the physical controller that requires practice. An example of when this open loop control was used is for the motion planning section. After giving the snakeboard an initial velocity, various changes were made using the GUI to determine how to slow it down and what affected the rate of deceleration.

– 42 –

The open loop periodic control simulation works by taking the controlled parameters (determined by the values in the GUI) and solving the equations of motion for a very small timestep. This is repeated continuously so that it appears to be in real time. Another option available to the user was enabling a “crumb trail”. At fixed intervals, a marker is left at the current position of the centre of mass of the snakeboard. This proved to be an effective way of tracking the path of the snakeboard and showing it to the user. An example of the output displayed to the user when the “crumb trail” is turned on is shown in Figure 64.

Figure 64 – “Crumb trail”

Although this simulation, using open loop periodic controls, was not used for the research done in this paper, it contributed indirectly. It proved invaluable as a tool to become familiar with the snakeboard as an abstract system and bridged the gap of the author‟s intuitive knowledge of snakeboards and nonholomonic systems. It is a good aid to visualising how the controlled parameters cause motion in the snakeboard and was inspiration for further work carried out in this paper.

– 43 –

9.

Bifurcation study

A bifurcation is said to occur when a small, gradual change applied to the input parameters causes an abrupt „qualitative‟ change in the behaviour of the system. For the case of the snakeboard system, this means a small alteration of the controlled parameters causes a qualitative change in the trajectory. Bifurcation theory was used to analyse the motion of the snakeboard in resonance gait to investigate its sensitivity and to gain a better understanding of its motion. The investigation was split between phase shifts and angular frequency changes. A MATLAB script was created that automated the process of trying each parameter combination. The trajectory of the snakeboard was plotted and saved as an image file to examine later. Each simulation in this section was run for 100 seconds. This relatively long duration was chosen to ensure the full behaviour of the snakeboard was captured rather than a small snapshot.

9.1.

Phase shift

9.1.1.

Rotor phase shift

The parameters in Table 27 were kept constant as 𝛽𝑟 was increased from 0 to 2 in increments of 0.01 radians. Parameter

Value

𝑎𝑓 𝑎𝑏 𝑎𝑟 𝜔𝑓 𝜔𝑏 𝜔𝑟 𝛽𝑓 𝛽𝑏

0.3 rad -0.3 rad 0.7 rad 1 rad/s 1 rad/s 1 rad/s 0 rad 0 rad

Table 27 – Constant parameters for rotor phase shift

A bifurcation point was found when 𝛽𝑟 = 1.5 radians which corresponds to the rotor being approximately 90° out of phase with the platforms. The trajectory is shown in Figure 65. The snakeboard travels a short distance in

– 44 –

the positive 𝑥 and 𝑦 direction before turning and travelling in the negative 𝑥 and 𝑦 direction.

Figure 65 – Bifurcation point at 𝜷𝒓 = 𝟏. 𝟓𝟖 radians

9.1.2.

Platform phase shift

The parameters in Table 28 were kept constant as 𝛽𝑓 was increased from 0 to 2 in increments of 0.01 radians. Another test was performed were the parameters in Table 29 were kept constant as 𝛽𝑏 was increased from 0 to 2 in increments of 0.01 radians. The purpose of running similar tests for 𝛽𝑓 and 𝛽𝑏 was to investigate whether or not there is a symmetrical property of the snakeboard where it does not matter if a phase change is made in the front or back platform. Parameter

Value

Parameter

Value

𝑎𝑓 𝑎𝑏 𝑎𝑟 𝜔𝑓 𝜔𝑏 𝜔𝑟 𝛽𝑏 𝛽𝑟

0.3 rad -0.3 rad 0.7 rad 1 rad/s 1 rad/s 1 rad/s 0 rad 0 rad

𝑎𝑓 𝑎𝑏 𝑎𝑟 𝜔𝑓 𝜔𝑏 𝜔𝑟 𝛽𝑓 𝛽𝑟

0.3 rad -0.3 rad 0.7 rad 1 rad/s 1 rad/s 1 rad/s 0 rad 0 rad

Table 28 – Constant parameters for first platform phase shift test

Table 29 – Constant parameters for second platform shift test

Although no bifurcation points were found, this test made is clear that it does matter which platform has the phase change applied to it. For example, a phase change of 0.69 is applied to the front platform (Figure 66) and the back platform (Figure 67).

– 45 –

Figure 66 – Trajectory when 𝜷𝒇 = 𝟎. 𝟎𝟔𝟗 rad

Figure 67 – Trajectory when 𝜷𝒃 = 𝟎. 𝟎𝟔𝟗 rad

The platform that is leading is more sensitive to changes than the trailing platform. When phase changes were applied to the back platform, the trajectory barely changed. However, phase changes applied to the front platform caused erratic behaviour and the snakeboard covered only half the distance.

9.2.

Angular Frequency

9.2.1.

Rotor angular frequency

The parameters in Table 30 were kept constant as 𝜔𝑟 was increased from 0 to 2 in increments of 0.01 rad/s. Parameter

Value

𝑎𝑓 𝑎𝑏 𝑎𝑟 𝜔𝑓 𝜔𝑏 𝛽𝑓 𝛽𝑏 𝛽𝑟

0.3 rad -0.3 rad 0.7 rad 1 rad/s 1 rad/s 0 rad 0 rad 0 rad

Table 30 – Constant parameters for rotor angular frequency test

– 46 –

It was found that although the trajectories varied greatly when the angular frequency way changed, there was no specific bifurcation point. 9.2.2.

Platform angular frequency

The parameters in Table 31 were kept constant as 𝜔𝑓 was increased from 0 to 2 in increments of 0.01 rad/s. Another test was performed where the parameters in Table 32 were kept constant as 𝜔𝑏 was increased from 0 to 2 in increments of 0.01 for the same reasons given in “Platform phase shift”. Parameter

Value

Parameter

Value

𝑎𝑓 𝑎𝑏 𝑎𝑟 𝜔𝑏 𝜔𝑟 𝛽𝑓 𝛽𝑏 𝛽𝑟

0.3 rad -0.3 rad 0.7 rad 1 rad/s 1 rad/s 0 rad 0 rad 0 rad

𝑎𝑓 𝑎𝑏 𝑎𝑟 𝜔𝑓 𝜔𝑟 𝛽𝑓 𝛽𝑏 𝛽𝑟

0.3 rad -0.3 rad 0.7 rad 1 rad/s 1 rad/s 0 rad 0 rad 0 rad

Table 31 – Constant parameters for first platform angular frequency test

Table 32 – Constant parameters for second platform angular frequency test

A bifurcation point was found for both the front and back platform at an angular frequency of 0.01 rad/s (Figure 70 and Figure 71). When either the front or back platform was stationary, the trajectory strongly resembled the resonance gait (Figure 68 and Figure 69). The implication of this is that motion can be achieved with only one rotating platform. However, the platform that is not rotating must stay perfectly still or the trajectory will be irregular.

Figure 68 – Trajectory when 𝝎𝒇 = 𝟎. 𝟎𝟏 rad/s

– 47 –

Figure 69 – Trajectory when 𝝎𝒃 = 𝟎. 𝟎𝟏 rad/s

Figure 70 – Trajectory when 𝝎𝒇 = 𝟎 rad/s

Figure 71 – Trajectory when 𝝎𝒃 = 𝟎 rad/s

Unlike the situation of front and back platforms being out of phase, there was no significant difference as to which platform was altered. An example of this is shown in the figures below. Figure 73 is the trajectory when 𝜔𝑓 = 1.21 rad/s and Figure 78 is the trajectory when 𝜔𝑏 = 1.21 rad/s. While not identical, they are very similar and this is representative of all the values of angular frequencies plotted.

Figure 72 – Trajectory when 𝝎𝒇 = 𝟏. 𝟐𝟏 rad/s

Figure 73 – Trajectory when 𝝎𝒃 = 𝟏. 𝟐𝟏 rad/s

– 48 –

10.

Motion Planning

Motion planning for nonholomonic systems is an area receiving a lot of academic attention. For this project, it was decided a basic motion planner should be created. The goal was to create an algorithm that would output controlled parameters that could enable the snakeboard to travel a path with a specified displacement in the 𝑥-axis in a given time. The resonance gait is used give a one-dimensional displacement and the snakeboard should start and end at rest. Since a resonance gait is used the angular frequencies for the platforms and rotors are assumed to equal (𝜔 = 𝜔𝑓 = 𝜔𝑏 = 𝜔𝑟 ) for motion planning. Two preliminary tests (described in “Speeding up and slowing down” and “Parameter to be changed”) were carried out to gain a better understanding of the snakeboard‟s behaviour in the process used for motion planning.

10.1.

Speeding up and slowing down

To prepare for motion planning tests were performed where the snakeboard was accelerated and then decelerated. These investigated how various distances could be travelled at different speeds and bring the snakeboard back to rest. 10.1.1. Test 1 The first test used the parameters in Table 33 for 30 seconds before changing to the parameters in Table 34 for the final 30 seconds. The trajectory is plotted in Figure 74 and the velocity against time in Figure 75. Parameter

Value

Parameter

Value

𝑎𝑓 𝑎𝑏 𝑎𝑟 𝜔 𝛽𝑓 𝛽𝑏 𝛽𝑟

0.3 rad -0.3 rad 0.7 rad 2 rad/s 0 rad 0 rad 0 rad

𝑎𝑓 𝑎𝑏 𝑎𝑟 𝜔 𝛽𝑓 𝛽𝑏 𝛽𝑟

0.3 rad -0.3 rad -0.7 rad 2 rad/s 0 rad 0 rad 0 rad

Table 33 – Parameters for initial period

Table 34 – Parameters for final period

– 49 –

Figure 74 – Trajectory during Test 1

Figure 75 – Velocity during Test 1

10.1.2. Test 2 The second test used the same parameters as Test 1 for the initial and final period (Table 33 and Table 34) but for 60 seconds each period. The trajectory is plotted in Figure 76 and the velocity against time in Figure 77

Figure 76 – Trajectory during Test 2

Figure 77 – Velocity during Test 2

– 50 –

The final velocity is 0.0135 m/s. This is further away from the ideal zero velocity than the previous test. From this, it was concluded that the longer the simulation was run for, the further the final velocity would be from zero. However, the final velocity is very low and could be considered negligible and is probably due to small errors in the simulation. For motion planning, either the time for deceleration could be increased so the snakeboard comes to completely rest or allow for a small error. In theory the final position of the snakeboard should have the same displacement in the 𝑦-axis as the initial position because of the symmetrical path created from accelerating and decelerating for the same length of time. Due to errors it is possible the displacement in the 𝑦-direction might not be exactly the same. However, the difference will be negligible compared to the displacement in the 𝑥-direction so it will be ignored. It was decided that in order to simplify the problem a tolerance of +/- 0.1 m/s would be applied.

10.2.

Parameter to be changed

From the resonant section there are three parameters that can be adjusted that would be suitable for varying the time taken in the motion planning algorithm by changing the acceleration of the snakeboard: rotor amplitude, platform amplitudes and the angular velocities. To test these parameters in the context of motion planning a test was performed for each parameter. For each increment the time taken to get to 0.5 m/s from rest and the distance travelled in the 𝑥-direction was recorded. The results from the tests in the previous section showed that if the same parameters were applied for the same length of time (except the rotor amplitude that was opposite) the snakeboard would come to rest. 10.2.1. Rotor Amplitude Test The parameters in Table 35 were kept constant as the rotor amplitude, 𝑎𝑟 , was increased from 0.2 radians to 1.5 radians in increments of 0.1 radians. The results are shown in Figure 78. – 51 –

Parameter

Value

𝑎𝑓 𝑎𝑏 𝜔 𝛽𝑓 𝛽𝑏 𝛽𝑟

0.3 rad -0.3 rad 1 rad/s 0 rad 0 rad 0 rad

Table 35 – Constant parameters for rotor amplitude test

Figure 78 – Time taken to reach 0.5 m/s for a range of rotor amplitudes

From Figure 78 it can be seen that there is an exponential relationship for the time taken and distance travelled as the rotor amplitude is increased. In this configuration an extremely large rotor amplitude would be needed to get to a velocity of 0.5 m/s in a time significantly less than 50 seconds. Likewise, an extremely large rotor amplitude would be needed to get to 0.5 m/s in a distance less than 10m. 10.2.2. Platform Amplitudes Test The parameters in Table 36 were kept constant as the platform amplitudes were increased from 0.2 to 1 radian. The results are shown in Figure 79. Parameter

Value

𝑎𝑟 𝜔 𝛽𝑓 𝛽𝑏 𝛽𝑟

0.7 rad 1 rad/s 0 rad 0 rad 0 rad

Table 36 – Constant parameters for platform amplitudes test

– 52 –

Figure 79 – Time taken to reach 0.5 m/s for a range of platofrm amplitudes

When the platform amplitude is varied from 0.2 to 0.6 radians, the relationship is similar to the rotor amplitude test. However, as was discussed in the resonance section the trajectory becomes unstable much faster when large platform amplitudes are used. 10.2.3. Velocity Test The parameters in Table 37 were kept constant as the angular velocities 𝜔 were increased from 0.8 to 2.0 rad/s in increments of 0.1 rad/s. The results are shown in Figure 80. Parameter

Value

𝑎𝑓 𝑎𝑏 𝑎𝑟 𝛽𝑓 𝛽𝑏 𝛽𝑟

0.3 rad -0.3 rad 0.7 rad 0 rad 0 rad 0 rad

Table 37 – Constant parameters for velocity test

– 53 –

Figure 80 – Time taken to reach 0.5 m/s for a range of angular frequencies

Figure 80 shows the time taken to reach 0.5 m/s. The distance travelled in the 𝑥-direction decreases exponentially as the angular velocity of the platforms and the rotor is increased. 10.2.4. Conclusions Varying the platform amplitudes can be immediately discounted for the purposes of this motion planning algorithm because the range they can be varied is very limited without causing instability. Whilst varying the rotor amplitude is suitable for slow accelerations, if a fast accelerating were required an extremely large amplitude would be necessary. In comparison, changing the angular velocities offers a large scope for motion planning. Therefore, it was decided that changing the angular velocity of the platforms and rotors is the best parameter to change.

– 54 –

10.3.

Algorithms

From the preliminary tests described previously, it was possible to begin developing a motion planning algorithm. All algorithms developed for this project are numerical solutions scripted in MATLAB. It was not originally intended to produce more than one algorithm. However, when the solutions were taking a long time to generate it was clear improvements needed to be made and the code optimised. A standard test was created (Table 38) to allow a comparison of the algorithms using the fixed control parameters (Table 39).

Requirement

Value

Displacement 30m Duration 100 seconds

Table 38 – Requirements to test motion planning algorithm

Parameter

Value

𝑎𝑓 𝑎𝑏 𝑎𝑟 𝛽𝑓 𝛽𝑏 𝛽𝑟

0.3 rad 0.3 rad -0.7 rad 0 rad 0 rad 0 rad

Table 39 – Fixed parameters

The time taken to produce a result at the given accuracy was measured. It is acknowledged that computing power available will determine this but since only a rough comparison was needed it would be suitable providing it was done on the same computer. All the algorithms in this paper have the snakeboard accelerating for half of the simulation time and decelerating for the second half.

– 55 –

10.3.1.

Version 1

The flowchart shown in Figure 81 describes the logic used in this algorithm. The algorithm begins by calculating the distance travelled at a minimal time period and ω. If the distance travelled was not sufficient, the period was increased until it was enough. If the current period was greater than the duration required, then the snakeboard will need to accelerate and decelerate faster to travel further. To do this 𝜔 is increased for the reasons given previously. Eventually an 𝜔 value will be found that allows the snakeboard to travel the required distance in the duration specified.

Get the xdisplacement, duration from the user. The current period is set to a minimum and 𝜔 to 0.1

Run the simulation for the current period at the current 𝜔 value

Is the calculated displacement greater than the required displacement?

Increase the period of the simulation by an increment of time

NO

YES

Increase 𝜔 by 𝜔increment. Reset current period to minimum

NO

Is the current period the same as the required time?

YES

Parameters successfully found, output to user

Figure 81 – Motion Planning (Version 1) Flowchart

– 56 –

Using the standard test, it was found that 𝜔 = 1.62 rad/s for 98 seconds caused the snakeboard to travel 15.4m. The time taken to reach the solution was 281 seconds. It was decided that the calculation time was unacceptably high, indicating the algorithm was very inefficient. The reason for this is that many time periods are evaluated for each value of 𝜔, even when that value is not a solution. For low values of 𝜔 this is a fundamental problem because it can require a very large period for the snakeboard to travel the required distance. Another problem with this algorithm is that it cannot begin at 𝜔 = 0 rad/s because no matter how long the simulation ran for it would not move. This means 𝜔 must start at an arbitrarily small value. If this happens to be larger than the solution then the results will be inaccurate. The accuracy of the solution can be improved by decreasing the time increment and the 𝜔 increment but this will increase the time taken to reach a solution.

– 57 –

10.3.2. Version 2 This version of the motion planning algorithm improved significantly on the previous one because it only evaluates each 𝜔 value once. Since the duration is specified, the simulation is evaluated for that period. The flowchart in Figure 82 describes the logic used in this algorithm. It begins with 𝜔 at zero and steadily increases 𝜔 until the snakeboard travels far enough. Each simulation is run for the specified duration.

Get the xdisplacement, duration from the user. Set 𝜔 to zero

Run the simulation for the specified duration at the current 𝜔 value

Is the calculated displacement greater than the required displacement?

NO

Increase 𝜔 by 𝜔increment

YES

Parameters successfully found, output to user Figure 82 – Motion Planning (Version 2) Flowchart

Using the standard test and 𝜔-increment = 0.01, it was found that 𝜔 = 1.57 for 100 seconds caused the snakeboard to travel 15.02m. The time taken to

– 58 –

reach the solution was 4.54 seconds. This algorithm is simpler than the previous one and offers a significant improvement because the solution time has been reduced by a significant factor, The accuracy of the result depends on the 𝜔 increment. A smaller increment will give increased accuracy at the expense of extra computational time because more iterations are needed. The standard test was run again but this time with an 𝜔-increment of 0.001. It was found that 𝜔 = 1.569 for 100 seconds caused the snakeboard to travel 15.0042m. The distance travelled is now closer to what was required, however, the computational time increased to 42.78 seconds. 10.3.3. Version 3 The previous version has a compromise between accuracy and speed. Because it is desirable to obtain an accurate solution in the least amount of time, a new algorithm was developed. The flowchart in Figure 83 describes the logic used in this algorithm. It is an extension of the previous version because it begins large 𝜔 increments until an estimate of the controlled parameters is found. It then refines 𝜔 with decreasing 𝜔-increments until a solution of sufficient accuracy is obtained. During testing, it was sometimes found that the solution would not converge to the accuracy required. This is probably due to small errors in the mathematical model and the non-linear nature of the system. When this happened, 𝜔-increment decreased exponentially towards zero. To prevent this from happening a minimum 𝜔-increment was defined and the algorithm aborted when this was reached.

– 59 –

Get the xdisplacement, duration and accuracy from the user. Set 𝜔 to zero

Run the simulation for the specified duration at the current 𝜔 value

Is the calculated displacement greater than the required displacement?

NO

Increase 𝜔 by 𝜔increment

YES

Parameters successfully found, output to user

NO

Go back to the previous 𝜔 value and reduce the 𝜔increment size

Did it overshoot the desired accuracy?

YES

Solution is not converging

YES

Is the 𝜔increment too small?

NO

Figure 83 – Motion Planning (Version 3) Flowchart

Using the standard test, it was found that 𝜔 = 1.556035041809082 rad/s for caused the snakeboard to travel 30.000004 m. The time taken to reach the solution was 2.63 seconds. The level of accuracy for 𝜔 needs to be extremely high to travel precise displacements. Again, the solution time has been reduced compared to the previous algorithms. Figure 85 shows how the iterative process efficiently converges upon the solution. As intended, the, 𝜔-increment decreases as the algorithm converges on the solution. – 60 –

Figure 84 – Iterative process

10.4.

Figure 85 – Output to user

Improvements

Only straight line motion planning is implemented in this algorithm. Although this is useful, full motion planning would require ability to orientate the board. This could be achieved by continuing the work done in the section titled “Rotation Gait”. The solutions generated by this motion planning algorithm have the snakeboard accelerate and decelerate for the same length of time. A better algorithm would have the snakeboard accelerate rapidly at the start with large omega and amplitudes before settling down to a constant velocity in a way similar to that of the results of the controller device. The stable resonant gait method could be employed to limit the amplitude of the oscillations about the 𝑥-axis. This would allow greater distances and eliminate the problems associated with instability. This algorithm does not take into account the physical limitations of a real world implementation, a heavy rotor can only accelerate so quickly. Ideally, the algorithm would compromise between different adjustments to the control parameters instead of only optimising 𝜔. For instance when the maximum amplitude of rotation for the platforms was reached the frequency would increase.

– 61 –

11.

Stable Resonance Gait

As mentioned previously, amplitudes in the 𝑦-direction can become large enough to cause the board to turn around in the resonant gait. This is particularly true when long simulation run times are used with large platform amplitudes. The purpose of the work in this section is create a script that covers the greatest distance in the 𝑥-distance in a given amount of time. It should be capable of long run times that are normally problematic when using the resonance gait. For the initial algorithm, the parameters that could change over time were 𝑎𝑓 , 𝑎𝑏 and 𝑎𝑟 . The platform angles were assumed equal and opposite (𝑎𝑝 = 𝑎𝑓 = −𝑎𝑏 ). Maximum and minimum values for 𝑎𝑝 and 𝑎𝑟 were assigned as parameters that could be controlled by the user. Because the resonance gait was used, the controlled parameters had the same frequency. As described in Figure 86, the controlled parameters can be split into “cycles”. Like wavelengths, their period depends on the frequency 𝜔 (for the purposes of this section, 𝜔 = 𝜔𝑓 = 𝜔𝑏 = 𝜔𝑟 ).

Cycle 1

Cycle 2

Figure 86 – “Cycles” in the controlled parameters for resonant gait

The algorithm works by evaluating every possible combination of amplitudes and determining which arrangement produces the greatest displacement in – 62 –

the 𝑥-direction for that cycle. The best combination of amplitudes is then run and the end conditions for that cycle are used for the beginning of the next cycle. This continues until the total duration of all the cycles is greater or equal to the given time. The accuracy number of combinations depended on the resolution in the range of 𝑎𝑝 and 𝑎𝑟 that were used. A small resolution produces results that are more accurate but takes longer to compute. The algorithm was tested using a range for 𝑎𝑝 from 0 to 0.7 in increments of 0.001 radians and a range of 𝑎𝑟 from 0 to 0.7 in increments of 0.001 radians. The angular frequency 𝜔 was fixed at 1 rad/s and the run time was 12𝜋 seconds.

Figure 87 – Trajectory during stable resonance gait test

Figure 88 – Controlled parameters during stable resonance gait test

– 63 –

Figure 89 – Platform and rotor amplitudes during stable resonance gait test

Figure 90 – Snakeboard velocity during stable resonance gait

The trajectory of the snakeboard (Figure 87) is not perpendicular to the 𝑥axis. It was decided that this was not a significant problem because an initial, small correction could be applied and the snakeboard‟s trajectory would be “straight”. The controlled parameters and 𝑎𝑝 and 𝑎𝑟 as functions of time are plotted in Figure 88 and Figure 89 respectively. When the snakeboard is initially at rest, the optimum platform amplitudes is 0.5 radians. This is similar to the amplitude that what was found in the section regarding the amplitude of the platforms in resonance gait. After a relatively optimum large platform angle in the first cycle the later cycles are much smaller at 0.02 and 0.029. This closely matches the author‟s experience in riding snakeboards and observing others. The best way to gain momentum is to begin with a large

– 64 –

rotation of the body (rotor) combined with large platform rotations. This gets the snakeboard moving before smaller oscillations are used. An initial spike in velocity (Figure 90) occurs before settling down. This simulation confirms that a big, initial “swing” is the optimum way of getting speed before settling down. A maximum rotor amplitude is always best for big 𝑥-displacement. This knowledge could be of use when designing robots that produce motion in a similar way to snakeboards.

– 65 –

12.

Conclusions 12.1.

Summary

By developing a mathematical model in MATLAB using the equations of motion for a snakeboard the author has been able to demonstrate some of it characteristics. The three basic gaits have been analysed to show how the direction of travel can be controlled. For each gait, the controlled parameters have been investigated to show how the behaviour can be altered. Using the knowledge gained from the basic gaits a motion planning algorithm was developed. It was found that accurate displacements required extremely precise control parameters that would be almost impossible for robots in the real world to replicate. To achieve reasonable accuracy a closed loop control mechanism would need to be implemented. Research in the bifurcation section showed that generally, the snakeboard does not exhibit many bifurcation points, instead, small parameter changes usually lead to small changes in the trajectory. An algorithm was created to resolve a problem highlighted in the resonance section that the trajectory became stable over time. This algorithm replicated the behaviour of a human rider on a snakeboard. Finally two methods were created to allow a user to control the simulation: open loop controls and a physical controller. The physical controller proved to be an interesting way of interacting with the simulation and was useful for testing different ideas.

– 66 –

12.2.

Evaluation

All the work done in this paper was based on a single MATLAB model. It would have been beneficial for real snakeboard behaviour to have been analysed. This would have given additional validation to the results. Efforts were made to ensure the code written was versatile and could be adapted for different aspects of the analysis. For this paper, it was soon realised a set of basic scripts should be created to plot the trajectory, controlled parameters and velocity for each simulation. This worked well because it ensured consistency throughout. For the author to continue this investigation a better understanding of the mathematics behind nonholomonic systems would be beneficial to enable a better understanding of the academic papers available.

12.3.

Further Work

To the best of the author‟s knowledge, nobody has yet to use motion capture technology on a human skateboarder. This could yield useful results for the field of robotics from analysing the way a human controls the board. Building a robotic snakeboard would be the next step to evaluating the motion planning algorithm developed in this paper. In addition, the algorithm could be extended beyond the linear point developed in this paper. This could be achieved using the rotation gait combined with the resonance gait. Finally, the stable resonance gait could be combined with the motion planning gait which would improve the latter greatly. The stable resonance gait provides a means for the snakeboard to achieve a reasonable speed quickly, rather than a slow acceleration.

– 67 –

13.

References

1. Nonholomonic mechanics and locomotion: the Snakeboard example. Lewis, Andrew D., et al. s.l. : IEEE International Conference on Robotics and Automation, 1994. 2. Nonholomonic Motion Planning: Steering Using Sinusoids. Murray, Richard. M. 5, s.l. : IEEE Transactions on Automatic Control, 1993, Vol. 38. 3. Minimum Control-Switch Motions for the Snakeboard: A Case Study in Kinematically Controllable Underactuated Systems. Ianniti, Stefano and Kevin, Lynch M. 4, s.l. : IEEE Transactions on Robotics, 2004, Vol. 20. 4. Kinematic controllability and motion planning for the snakeboard. Bullo, F and Lewis, A D. s.l. : IEEE Transaction on Robotics, 2003, Vol. 19. 5. Some New Robust Pseudo Forward and Rotation Gaits for the Snakeboard. Asnafi, A. and Mahzoon, M. 5, s.l. : Scientia Iranica, Vol. 15. 6. Veltman, Eddy. Design and realization of an experimental snakeboard. s.l. : University of Twente, 2004. 7. Further Development of the Mathematical Model of a Snakeboard. Kuleshow, A. S. Moscow : s.n., 2007.

– 68 –

Appendix I – Timestep

Appendix II – MATLAB Code 1. Differential Equation Solver % % % % % %

Title: Snakeboard Differential Equation Solver Filename: snakeboard_differential_equations.m Author: Jonathan Jamieson Date: September 2011 – April 2012 Purpose: This function solves the equations of motion that govern a snakeboard's motion over a specified timespan

function snakeboard=snakeboard_differential_equations(t,x,m,l,J,Jr,Jp,a1,a2,ar,omega_1,omega_2 ,omega_r,Bf,Bb,Br) k_squared = ((J+Jr+2*Jp)/(m*l^2)); d1 = (Jr/m*l); d2 = (Jp/m*l); psi_1 = (a1*sin(omega_1*t+Bf)-a2*sin(omega_2*t+Bb)); psi_1_dot = (omega_1*a1*cos(omega_1*t)-omega_2*a2*cos(omega_2*t+Bb)); psi_1_2dot = (((omega_1^2)*omega_1*sin(omega_1*t))+((omega_2^2)*omega_2*sin(omega_2*t+Bb))); psi_2 = (a1*sin(omega_1*t+Bf)+a2*sin(omega_2*t+Bb)); psi_2_dot =(omega_1*a1*cos(omega_1*t+Bf)+omega_2*a2*cos(omega_2*t+Bb)); psi_2_2dot = ((-(omega_1^2)*omega_1*sin(omega_1*t+Bf))((omega_2^2)*omega_2*sin(omega_2*t+Bb))); delta = (ar*sin(omega_r*t+Br)); delta_dot = (omega_r*ar*cos(omega_r*t+Br)); delta_2dot = (-((omega_r)^2)*ar*sin(omega_r*t+Br)); snakeboard=[((x(4)*cos(x(3)))(x(4)*(sin((psi_2)))*sin(x(3)))/(cos((psi_1))+cos((psi_2)))); ((x(4)*sin(x(3)))+(x(4)*(sin((psi_2)))*cos(x(3)))/(cos((psi_1))+cos((psi_2)))); ((x(4)*sin((psi_1)))/(l*(cos((psi_1))+cos((psi_2))))) (((((d1*((delta_2dot))+d2*((psi_2_2dot)))*sin((psi_1)))/(cos((psi_1))+cos((psi_2))))(x(4)*((((psi_2_dot)*sin((psi_2))*cos((psi_2))+k_squared*(psi_1_dot)*sin((psi_1))*cos ((psi_1)))/(cos((psi_1))+cos((psi_2)))^2)+((((psi_1_dot)*sin((psi_1)))+((psi_2_dot)*s in((psi_2))))/((cos((psi_1))+cos((psi_2)))^3))*(((sin((psi_2)))^2)+(k_squared*(sin((p si_1)))^2)))))/(1+((((psi_2))^2)+(k_squared)*((psi_1))^2)/((cos((psi_1))+cos((psi_2)) )^2)))];

2. General Snakeboard Motion % % % % %

Title: General Snakeboard Motion Filename: general.m Author: Jonathan Jamieson Date: September 2011 – April 2012 Purpose: Change parameters to experiment with various gaits

%% Clear Previous figures, variables etc clc;

clear all; close all; %% a1 a2 ar

Control Parameters = +0.3; = -0.3; = 0.7;

omega_1 = 1; omega_2 = 1; omega_r = 1; Bf = 0; Bb = 0; Br = 0; %% Define initial conditions and Timespan of Calculation x0=[0;0;0;0]; % Define initial conditions (x position, y position, theta, velocity); initial_time = 0; % Starting time final_time = 500; % Finishing time timestep = 0.1; % The resolution of the differentiation %% Define snakeboard properties m = 75; % kg, mass of snakeboard and rider l = 0.285; % m, half the length of the snakeboard J = 0.22; % moment of inertia of crossbar and platforms Jr = 14; % moment of inertia of rotor Jp = 0.013; % moment of inertia of platforms %% Simulation Properties tspan=initial_time:timestep:final_time; %Define timespan (period to integrate)

%% Solve coupled differential equations (equations of motion of the system) [t,x]=ode45(@snakeboard_differential_equations,tspan,x0,[],m,l,J,Jr,Jp,a1,a2,ar,omega _1,omega_2,omega_r,Bf,Bb,Br);

%% Calculate angles over time span phi_f = a1*sin(omega_1*tspan+Bf); % Front Platorm Angle phi_b = a2*sin(omega_2*tspan+Bb); % Back Platform Angle delta = ar*sin(omega_r*tspan+Br); % Rotor Angle

%% Extract values from solver xpos = x(:,1); ypos = x(:,2); theta = x(:,3); v= x(:,4); %% Write state and position of the snakeboard as a function of time to a matrix snakeboard = zeros(numel(tspan),8); snakeboard(:,1) snakeboard(:,2) snakeboard(:,3) snakeboard(:,4) snakeboard(:,5) snakeboard(:,6) snakeboard(:,7) snakeboard(:,8)

= = = = = = = =

tspan; xpos; ypos; theta; v; phi_f; phi_b; delta;

3. Time Step Test % % % % %

Title: Time Step Test Filename: time_step_test.m Author: Jonathan Jamieson Date: September 2011 – April 2012 Purpose: Used to test the behaviour of the ODE45 solver

%% Clear Previous figures, variables etc clc; clear all; close all; a1 = +0.3; a2 = -0.3; ar = 0.3; omega_1 = 1; omega_2 = 1; omega_r = 1; Bf = 0; Bb = 0; Br = 0; %% Define initial conditions and Timespan of Calculation x0=[0;0;0;0]; % Define initial conditions (x position, y position, theta, velocity); initial_time = 0; % Starting time final_time = 30; % Finishing time timestep = 5; % The resolution of the differentiation5 %% Define snakeboard properties m = 75; % kg, mass of snakeboard and rider l = 0.285; % m, half the length of the snakeboard J = 0.22; % kg.m^2, moment of inertia of crossbar and platforms Jr = 14; % moment of inertia of rotor Jp = 0.013; % moment of inertia of platforms scrsz = get(0,'ScreenSize'); %% Simulation Properties tspan=initial_time:timestep:final_time; %Define timespan (period to integrate)

%% Solve coupled differential equations (equations of motion of the system) tic; [t,x]=ode45(@snakeboard_differential_equations,tspan,x0,[],m,l,J,Jr,Jp,a1,a2,ar,omega _1,omega_2,omega_r,Bf,Bb,Br); ode_time = toc; %% Calculate angles over time span phi_f = a1*sin(omega_1*tspan+Bf); % Front Platorm Angle phi_b = a2*sin(omega_2*tspan+Bb); % Back Platform Angle delta = ar*sin(omega_r*tspan+Br); % Rotor Angle %% Control Variables over time psi_1 = (a1*sin(omega_1*t+Bf)-a2*sin(omega_2*t+Bb)); psi_2 = (a1*sin(omega_1*t+Bf)+a2*sin(omega_2*t+Bb)); delta_control = (ar*sin(omega_r*t+Br));

%% Create Plots

f = figure('Position',[scrsz(1)+scrsz(3)/4 scrsz(2)+scrsz(4)/6 scrsz(3)/2 scrsz(4)/1.5]); name = ['Trajectory when time step = ',num2str(timestep)]; plot(x(:,1),x(:,2)),xlabel('x (m)'),ylabel('y (m)'),title(name); xpos = x(:,1); final_x = xpos(end)

4. Define Surface Plot % % % % % %

Title: Define Surface Plot Filename: surface.m Author: Jonathan Jamieson Date: September 2011 – April 2012 Purpose: Calls the surface plot function and gives it the input parameters

%% Clear Previous figures, variables etc clc; clear all; close all; %% Control Parameters a1 = 0.3; a2 = -0.3; ar = 0.7; omega_1 = 1; omega_2 = 1; omega_r = 1; Bf = 0; Bb = 0; Br = 0; %% Define initial conditions and Timespan of Calculation x0=[0;0;0;0]; % Define initial conditions (x position, y position, theta, velocity); initial_time = 0; % Starting time final_time = 200; % Finishing time timestep = 0.1; % The resolution of the differentiation m = 75; % mass of snakeboard and rider l = 0.285; % half the length of the snakeboard J = 0.22; % moment of inertia of crossbar and platforms Jr = 14; % moment of inertia of rotor Jp = 0.013; % moment of inertia of platforms X = 0:0.05:2; Y = 0:0.05:2; tspan=initial_time:timestep:final_time; %Define timespan (period to integrate) Z=surface_function(X,Y,m,l,J,Jr,Jp,a1,a2,ar,omega_1,omega_2,omega_r,Bf,Bb,Br,tspan,x0 ); % Calculate displacement for each combination surf(X,Y,Z), shading interp,xlabel('X Label'),ylabel('Y Label'),zlabel('Displacement in the x-direction (m)'); % Draw surface plot

5. Surface Plot Function % % % % % %

Title: Surface Plot Function Filename: surface_function.m Author: Jonathan Jamieson Date: September 2011 – April 2012 Purpose: This function creates a matrix (z) that consists of how far the snakeboard has gone in the x-direction

function z = surface_function(X,Y,m,l,J,Jr,Jp,a1,a2,ar,omega_1,omega_2,omega_r,Bf,Bb,Br,tspan,x0); z = zeros(numel(Y),numel(X)); % Predefine size of matrix

for i_x = 1:numel(X) for i_y = 1:numel(Y) % Can be angular frequencies, amplitude, or phase shift omega_1 = X(i_x); omega_2 = X(i_x); omega_r = Y(i_y); [t,x]=ode45(@snakeboard_differential_equations,tspan,x0,[],m,l,J,Jr,Jp,a1,a2,ar,omega _1,omega_2,omega_r,Bf,Bb,Br); xpos = x(:,1); x_distance = xpos(end); z(i_y,i_x) = x_distance; end end

6. Parameter variation & rate of change tests % % % % %

Title: Parameter varitation & rate of change tests Filename: parameter_variation.m Author: Jonathan Jamieson Date: September 2011 – April 2012 Purpose: Defines inputs for the function and plots the results

%% Clear Previous figures, variables etc clc; clear all; close all; %% a1 a2 ar

Control Parameters = 0.3; = -0.3; = 0.7;

omega_1 = 1; omega_2 = 1; omega_r = 1; Bf = 0; Bb = 0; Br = 0;

%% Define initial conditions and Timespan of Calculation x0=[0;0;0;0]; % Define initial conditions (x position, y position, theta, velocity); initial_time = 0; % Starting time final_time = 200; % Finishing time timestep = 0.1; % The resolution of the differentiation

m = 75; % mass of snakeboard and rider l = 0.285; % half the length of the snakeboard J = 0.22; % moment of inertia of crossbar and platforms Jr = 14; % moment of inertia of rotor Jp = 0.013; % moment of inertia of platforms X = 0:0.01:0.7; tspan=initial_time:timestep:final_time; % Define timespan (period to integrate) Z=calculate_measured_parameter(X,m,l,J,Jr,Jp,a1,a2,ar,omega_1,omega_2,omega_r,Bf,Bb,B r,tspan,x0); plot(Z(:,1),Z(:,2)),xlabel('X Label'),ylabel('Y Label ');

7. Calculate Measured Parameter % % % % %

Title: Calculated measured parameter Filename: calculate_measured_parameter.m Author: Jonathan Jamieson Date: September 2011 – April 2012 Purpose: Calculates results for each combination of the parameters

function z = calculate_measured_parameter(X,m,l,J,Jr,Jp,a1,a2,ar,omega_1,omega_2,omega_r,Bf,Bb,Br, tspan,x0); z = zeros(numel(X),2); % Predefine size of matrix for i_x = 1:numel(X) a1 = X(i_x); a2 = -X(i_x); [t,x]=ode45(@snakeboard_differential_equations,tspan,x0,[],m,l,J,Jr,Jp,a1,a2,ar,omega _1,omega_2,omega_r,Bf,Bb,Br); xpos = x(:,1); x_distance = xpos(end); z(i_x,2) = x_distance; z(i_x,1) = a1; end

8. Physical simulation controller % % % % % %

Title: Physical simulation controller Filename: controller.m Author: Jonathan Jamieson Date: September 2011 – April 2012 Purpose: This function solves the equations of motion that govern a snakeboard's motion over a specified timespan

%% Clear Previous figures, variables etc clc; clear all; close all; %% Define initial conditions and Timespan of Calculation x0=[0;0;0;0]; % Define initial conditions (x position, y position, theta, velocity); m = 75; % kg, mass of snakeboard and rider l = 0.285; % m, half the length of the snakeboard J = 0.22; % kg.m^2, moment of inertia of crossbar and platforms Jr = 14; % moment of inertia of rotor Jp = 0.013; % moment of inertia of platforms % Get Screen size scrsz = get(0,'ScreenSize'); % Serial Settings s = serial ('COM3'); s.BaudRate = 4800; s.Terminator = 'cr'; fopen (s); %% Simulation Properties figure; axis_size = 100; grid_lines = 1; ax = axes('XLim',[-axis_size axis_size],'YLim',[-axis_size axis_size],'ZLim',[-1.5 1.5]); view(2); grid on; [Xg,Yg] = meshgrid(-axis_size:grid_lines:axis_size, -axis_size:grid_lines:axis_size); Zg =.00001*Xg+.0001*Yg; surf(Xg,Yg,Zg) colormap white axis equal grid on; set(gcf,'Color','white') %% Define snakeboard graphical output %Everything is relative to half the length of snakeboard (l) % Crossbar [xb yb zb]=cylinder([l l]); g(1) = surface(1*xb,0.07*yb,0.07*zb+0.06,'FaceColor','black'); tb = hgtransform('Parent',ax); set(g,'Parent',tb) % Rotor [xr yr zr]=cylinder([l l]); h(1) = surface(0.5*xr,0.05*yr,0.07*zr+0.08,'FaceColor','green'); tr = hgtransform('Parent',ax);

set(h,'Parent',tr) % Front Platform [xp1 yp1 zp1]=cylinder([l l]); i(1) = surface(0.07*xp1,0.5*yp1,0.07*zp1,'FaceColor','yellow'); tp1 = hgtransform('Parent',ax); set(i,'Parent',tp1) % Back Platform [xp2 yp2 zp2]=cylinder([l l]); j(1) = surface(0.07*xp2,0.5*yp2,0.07*zp2,'FaceColor','blue'); tp2 = hgtransform('Parent',ax); set(j,'Parent',tp2) set(gcf,'Renderer','opengl') drawnow psi_1_old = 0; psi_2_old = 0; delta_old = 0; initial_time = 0; tic; snakeboard = zeros(2000,17); for t_run=1:1:2000, pot_values_raw = fscanf(s); pot_values_raw = fscanf(s); pot_values = str2num(pot_values_raw); phi_f = (0*pi/180)+((40*pi/180-0*pi/180)/(776-548))*(pot_values(1)-548); phi_b = (0*pi/180)+((40*pi/180-0*pi/180)/(785-557))*(pot_values(2)-557); delta = (0*pi/180)+((60*pi/180-0*pi/180)/(570-420))*(pot_values(3)-420); if phi_f > 80*pi/180 phi_f = 80*pi/180; elseif phi_f < -80*pi/180 phi_f = -80*pi/180; end if phi_b > 80*pi/180 phi_b = 80*pi/180; elseif phi_b < -80*pi/180 phi_b = -80*pi/180; end

if delta > 90*pi/180 delta = 90*pi/180; elseif delta < -90*pi/180 delta = -90*pi/180; end

time_change = toc; final_time = initial_time+time_change; psi_1 = phi_f - phi_b; psi_2 = phi_f + phi_b; psi_1_dot = (psi_1-psi_1_old)/time_change;

psi_1_2dot = psi_1_dot/time_change; psi_2_dot = (psi_2-psi_2_old)/time_change; psi_2_2dot = psi_2_dot/time_change; delta_dot = (delta-delta_old)/time_change; delta_2dot = delta_dot/time_change; tic tspan=initial_time:time_change/2:final_time; initial_time = final_time; psi_1_old = psi_1; psi_2_old = psi_2; delta_old = delta; %% Solve coupled differential equations (equations of motion of the system) [t,x]=ode45(@snakeboard_differential_equations_controller,tspan,x0,[],m,l,J,Jr,Jp,psi _1,psi_2,delta,psi_1_dot,psi_2_dot,delta_dot,psi_1_2dot,psi_2_2dot,delta_2dot); xpos = x(:,1); ypos = x(:,2); theta = x(:,3); velocity = x(:,4); snakeboard(t_run,:) = [initial_time xpos(end) ypos(end) theta(end) velocity(end) phi_f phi_b delta psi_1 psi_2 delta psi_1_dot psi_2_dot delta_dot psi_1_2dot psi_2_2dot delta_2dot]; %% Animate snakeboard

% Draw crossbar trans = makehgtform('translate',[xpos(end) ypos(end) 0]); rotz = makehgtform ('zrotate',theta(end)); set(tb,'Matrix',trans*rotz); % Draw Rotor trans = makehgtform('translate',[xpos(end) ypos(end) 0]); rotz = makehgtform ('zrotate',(theta(end)+delta)); set(tr,'Matrix',trans*rotz); % Draw Front Platform trans = makehgtform('translate',[(xpos(end)+l*cos(theta(end))) (ypos(end)+l*sin(theta(end))) 0]); rotz = makehgtform ('zrotate',(theta(end)+phi_f)); set(tp1,'Matrix',trans*rotz); % Draw Back Platfrom trans = makehgtform('translate',[(xpos(end)-l*cos(theta(end))) (ypos(end)l*sin(theta(end))) 0]); rotz = makehgtform ('zrotate',(theta(end)+phi_b)); set(tp2,'Matrix',trans*rotz); % Camera Position camzoom(1); camtarget([xpos(end),ypos(end),0]); campos([xpos(end)-0,ypos(end)+7,15]); set(gcf,'Renderer','opengl') drawnow %% Obtain results from previous differentiation and set that as the new starting value

x0 = x(end,:); x0 = x0.'; end datetime=datestr(now); datetime=strrep(datetime,':','_'); %Replace colon with underscore datetime=strrep(datetime,'-','_');%Replace minus sign with underscore datetime=strrep(datetime,' ','_');%Replace space with underscore save(datetime, 'snakeboard'); fclose (s) delete (s)

9. Differential equation solver for physical simulation controller % % % % % %

Title: Differential equation solver for physical simulation controller Filename:snakeboard_differential_equations_controller.m Author: Jonathan Jamieson Date: September 2011 – April 2012 Purpose: Solves the equations of motion for the simulation controller

function snakeboard=snakeboard_differential_equations_controller (t,x,m,l,J,Jr,Jp,psi_1,psi_2,delta,psi_1_dot,psi_2_dot,delta_dot,psi_1_2dot,psi_2_2do t,delta_2dot) k_squared = ((J+Jr+2*Jp)/(m*l^2)); d1 = (Jr/m*l); d2 = (Jp/m*l);

snakeboard=[((x(4)*cos(x(3)))(x(4)*(sin((psi_2)))*sin(x(3)))/(cos((psi_1))+cos((psi_2)))); ((x(4)*sin(x(3)))+(x(4)*(sin((psi_2)))*cos(x(3)))/(cos((psi_1))+cos((psi_2)))); ((x(4)*sin((psi_1)))/(l*(cos((psi_1))+cos((psi_2))))) (((((d1*((delta_2dot))+d2*((psi_2_2dot)))*sin((psi_1)))/(cos((psi_1))+cos((psi_2))))(x(4)*((((psi_2_dot)*sin((psi_2))*cos((psi_2))+k_squared*(psi_1_dot)*sin((psi_1))*cos ((psi_1)))/(cos((psi_1))+cos((psi_2)))^2)+((((psi_1_dot)*sin((psi_1)))+((psi_2_dot)*s in((psi_2))))/((cos((psi_1))+cos((psi_2)))^3))*(((sin((psi_2)))^2)+(k_squared*(sin((p si_1)))^2)))))/(1+((((psi_2))^2)+(k_squared)*((psi_1))^2)/((cos((psi_1))+cos((psi_2)) )^2)))];

10. % % % % %

Open loop periodic controls

Title: Open loop periodic controls Filename:open_loop.m Author: Jonathan Jamieson Date: September 2011 – April 2012 Purpose: Open loop control code

%% Clear Previous figures, variables etc clc; clear all; close all; %% Open GUI fig = openfig('gui_1.fig','reuse');

% reuse or new

%% Define initial conditions and Timespan of Calculation x=[0 0 0 0]; % Define initial conditions (x position, y position, theta, velocity); timestep = 0.25; % The resolution of the differentiation time_interval = 0.75; % Period of time to integrate over %% Define snakeboard properties m = 75; % mass of snakeboard and rider l = 0.285; % half the length of the snakeboard J = 0.22; % moment of inertia of crossbar and platforms Jr = 14; % moment of inertia of rotor Jp = 0.013; % moment of inertia of platforms %Animation Properties axis_size = 100; % How large to make the output grid grid_lines = 0.5; % Distance between grid lines crumb_frequency = 5; % How often "crumbs" are dropped. The smaller the number the more often they are drooped. Must be an integer.

%% Simulation Properties tspan = -timestep; % This ensures the simulation starts at 0 seconds animation_pause = 0.01; % Increase or reduce as necessary to give real time result

%% Create animation output window crumb_counter = 0; % Crumbs can be dropped every so often so the path of the snakeboard can be observed ax = axes('XLim',[-axis_size axis_size],'YLim',[-axis_size axis_size],'ZLim',[-1.5 1.5]); view(2); grid on;

[Xg,Yg] = meshgrid(-axis_size:grid_lines:axis_size, -axis_size:grid_lines:axis_size); Zg =.00001*Xg+.0001*Yg; surf(Xg,Yg,Zg) colormap white axis equal grid on; set(gcf,'Color','white') scrsz = get(0,'ScreenSize'); set(gcf,'Position',[scrsz(1)+scrsz(3)/8 scrsz(2)+scrsz(4)/8 scrsz(3)*3/4 scrsz(4)*3/4]); % This sets the size and the position of the animation output window

%% Define snakeboard graphical output %Everything is relative to half the length of snakeboard (l) % Crossbar [xb yb zb]=cylinder([l l]); g(1) = surface(1*xb,0.07*yb,0.07*zb+0.06,'FaceColor','black'); tb = hgtransform('Parent',ax); set(g,'Parent',tb) % Rotor [xr yr zr]=cylinder([l l]); h(1) = surface(0.5*xr,0.05*yr,0.07*zr+0.08,'FaceColor','green'); tr = hgtransform('Parent',ax); set(h,'Parent',tr)

% Front Platform [xp1 yp1 zp1]=cylinder([l l]); i(1) = surface(0.07*xp1,0.5*yp1,0.07*zp1,'FaceColor','yellow'); tp1 = hgtransform('Parent',ax); set(i,'Parent',tp1) % Back Platform [xp2 yp2 zp2]=cylinder([l l]); j(1) = surface(0.07*xp2,0.5*yp2,0.07*zp2,'FaceColor','blue'); tp2 = hgtransform('Parent',ax); set(j,'Parent',tp2) set(gcf,'Renderer','opengl') drawnow while(2 > 1)

tspan= tspan(end)+timestep:timestep:tspan(end)+time_interval; %Define timespan (period to integrate) %% Obtain results from previous differentiation and set that as the new starting value x0 = x(end,:); x0 = x0.'; %% Get GUI data slider_input = gui_1; slider_input = guidata(slider_input); a1 = get(slider_input.a1,'Value'); a2 = get(slider_input.a2,'Value'); ar = get(slider_input.ar,'Value'); omega_1 = get(slider_input.omega1,'Value'); omega_2 = get(slider_input.omega2,'Value'); omega_r = get(slider_input.omegar,'Value'); Bf = get(slider_input.Bf,'Value'); Bb = get(slider_input.Bb,'Value'); Br = get(slider_input.Br,'Value');

%% Solve coupled differential equations (equations of motion of the system) [t,x]=ode45(@snakeboard_differential_equations,tspan,x0,[],m,l,J,Jr,Jp,a1,a2,ar,omega _1,omega_2,omega_r,Bf,Bb,Br);

%% Calculate angles over time span phi_f = a1*sin(omega_1*tspan+Bf); % Front Platorm Angle phi_b = a2*sin(omega_2*tspan+Bb); % Back Platform Angle delta = ar*sin(omega_r*tspan+Br); % Rotor Angle %% Control Variables over time psi_1 = (a1*sin(omega_1*t+Bf)-a2*sin(omega_2*t+Bb)); psi_2 = (a1*sin(omega_1*t+Bf)+a2*sin(omega_2*t+Bb)); delta_control = (ar*sin(omega_r*t+Br));

%% Extract time series data for x position, y position, angle and velocity of snakeboard xpos = x(:,1); ypos = x(:,2);

theta = x(:,3); velocity = x(:,4);

[xc yc zc]=cylinder([.01 .01]); %% Animate snakeboard for i = 1:numel(xpos) % Draw crossbar trans = makehgtform('translate',[xpos(i) ypos(i) 0]); rotz = makehgtform ('zrotate',theta(i)); set(tb,'Matrix',trans*rotz); % Draw Rotor trans = makehgtform('translate',[xpos(i) ypos(i) 0]); rotz = makehgtform ('zrotate',(theta(i))+delta(i)); set(tr,'Matrix',trans*rotz); % Draw Front Platform trans = makehgtform('translate',[(xpos(i)+l*cos(theta(i))) (ypos(i)+l*sin(theta(i))) 0]); rotz = makehgtform ('zrotate',(theta(i))+phi_f(i)); set(tp1,'Matrix',trans*rotz); % Draw Back Platfrom trans = makehgtform('translate',[(xpos(i)-l*cos(theta(i))) (ypos(i)l*sin(theta(i))) 0]); rotz = makehgtform ('zrotate',((theta(i))+phi_b(i))); set(tp2,'Matrix',trans*rotz);

end % Camera Position camzoom(1); camtarget([xpos(i),ypos(i),0]); campos([xpos(i)-7,ypos(i)-7,10]);

% pause(animation_pause) % Crumb Trail crumb_counter = crumb_counter+1; if crumb_counter == crumb_frequency r(1) = surface(xc+xpos(i),yc+ypos(i),.05*zc,'FaceColor','black'); rr = hgtransform('Parent',ax); set(r,'Parent',rr) set(gcf,'Renderer','opengl') drawnow hold on crumb_counter = 0; end end

11. % % % % %

Bifurcation Test

Title: Bifurcation test Filename: bifurcation.m Author: Jonathan Jamieson Date: September 2011 – April 2012 Purpose: Performs the bifurcation tests

%% Clear Previous figures, variables etc clc; clear all; close all; %% Control Parameters a1 = 0.3; a2 = -0.3; ar = 0.7; omega_1 = 1; omega_2 = 1; omega_r = 1; Bf = 0; Bb = 0; Br = 0; %% Define initial conditions and Timespan of Calculation x0=[0;0;0;0]; % Define initial conditions (x position, y position, theta, velocity); initial_time = 0; % Starting time final_time = 100; % Finishing time timestep = 0.1; % The resolution of the differentiation %% Define snakeboard properties m = 75; % kg, mass of snakeboard and rider l = 0.285; % m, half the length of the snakeboard J = 0.22; % kg.m^2, moment of inertia of crossbar and platforms Jr = 14; % moment of inertia of rotor Jp = 0.013; % moment of inertia of platforms %% Simulation Properties tspan=initial_time:timestep:final_time; %Define timespan (period to integrate) for omega_2 = 0:0.01:2 %% Solve coupled differential equations (equations of motion of the system) [t,x]=ode45(@snakeboard_differential,tspan,x0,[],m,l,J,Jr,Jp,a1,a2,ar,omega_1,omega_2 ,omega_r,Bf,Bb,Br); %% Calculate angles over time span phi_f = a1*sin(omega_1*tspan+Bf); % Front Platorm Angle phi_b = a2*sin(omega_2*tspan+Bb); % Back Platform Angle delta = ar*sin(omega_r*tspan+Br); % Rotor Angle %% Control Variables over time psi_1 = (a1*sin(omega_1*t+Bf)-a2*sin(omega_2*t+Bb)); psi_2 = (a1*sin(omega_1*t+Bf)+a2*sin(omega_2*t+Bb)); delta_control = (ar*sin(omega_r*t+Br)); %% Create Plots scrsz = get(0,'ScreenSize'); f = figure('Position',[scrsz(1)+scrsz(3)/4 scrsz(2)+scrsz(4)/6 scrsz(3)/2 scrsz(4)/1.5]); subplot(2,1,1) plot(t,phi_f,t,phi_b,t,delta),xlabel('Time (seconds)'),ylabel('Angle (rads)'),legend('Front Platform Angle','Back Platform Angle','Rotor Angle',1),title('Controlled Parameters'); subplot(2,1,2) plot(x(:,1),x(:,2)),xlabel('x (m)'),ylabel('y (m)'),title('Trajectory');

str = ['omega_f=',num2str(omega_1),', omega_b=',num2str(omega_2),', omega_r=',num2str(omega_r),]; v = get(gca); lh = line([0 0 NaN v.XLim],[v.YLim NaN 0 0 ]); set(lh,'Color',[.25 .25 .25],'LineStyle',':'); uicontrol('Style', 'text',... 'String', str,... %replace something with the text you want 'Units','normalized',... 'Position', [0.8 0.65 0.1 0.1]); filename = sprintf('a1_%4.3f-a2_%4.3f-ar_%4.3f-omega1_%4.3f-omega2_%4.3fomegar_%4.3f-Bf_%4.3f-Bb_%4.3fBr_%4.3f.png',a1,a2,ar,omega_1,omega_2,omega_r,Bf,Bb,Br); print(f,'-dpng',filename); filename = sprintf('a1_%4.3f-a2_%4.3f-ar_%4.3f-omega1_%4.3f-omega2_%4.3fomegar_%4.3f-Bf_%4.3f-Bb_%4.3fBr_%4.3f.png',a1,a2,ar,omega_1,omega_2,omega_r,Bf,Bb,Br);

dpi = 200; mag = dpi / get(0, 'ScreenPixelsPerInch'); set(gca, 'Color', 'none'); % Sets axes background export_fig(gcf,'png','bounds',filename,'-transparent', sprintf('-m%g', mag)); close; end

12. % % % % %

Motion Planning (Version 3)

Title: Motion Planning Filename:motion_planning.m Author: Jonathan Jamieson Date: September 2011 – April 2012 Purpose: Calculates parameters to satisfy user’s requirements

%% Clear Previous figures, variables etc clc; clear all; close all;

%% Prompt user for input prompt = {'Enter x-displacement (metres):','Enter duration (seconds):','Enter accuracy required (m):'}; dlg_title = 'Snakeboard Motion Planning'; num_lines = 1; def = {'30','60','0.01'}; answer = inputdlg(prompt,dlg_title,num_lines,def); tic; x_displacement_required = str2num(answer{1}); duration = str2num(answer{2}); x_displacement_accuraccy = str2num(answer{3}); %% Control Parameters a1 = 0.3; a2 = -0.3; ar = 0.7; omega = 0; % Starting omega omega_increment = 0.5; minimum_omega_increment = 0.1*10^-10; monitor = [0 0 0];

iteration = 0;

Bf = 0; Bb = 0; Br = 0; achieved_distance = 0; x_distance = 0; %% Define initial conditions and Timespan of Calculation x0=[0;0;0;0]; % Define initial conditions (x position, y position, theta, velocity); time_step = 0.1; % The resolution of the differentiation m = 75; % kg, mass of snakeboard and rider l = 0.285; % m, half the length of the snakeboard J = 0.22; % kg.m^2, moment of inertia of crossbar and platforms Jr = 14; % moment of inertia of rotor Jp = 0.013; % moment of inertia of platforms tspan1 = 0:time_step:duration/2; tspan2 = duration/2:time_step:duration; while achieved_distance < 1 omega = omega + omega_increment; current_monitor = [iteration omega x_distance]; monitor = [monitor;current_monitor]; iteration = iteration + 1; x0=[0;0;0;0]; [t,x]=ode45(@snakeboard_differential_equations,tspan1,x0,[],m,l,J,Jr,Jp,a1,a2,ar,omeg a,omega,omega,Bf,Bb,Br); x0 = x(end,:); x0 = x0.'; xpos = x(:,1); ypos = x(:,2); theta = x(:,3); v = x(:,4); [t,x]=ode45(@snakeboard_differential_equations,tspan2,x0,[],m,l,J,Jr,Jp,a1,a2,ar,omega,omega,omega,Bf,Bb,Br); xpos = [xpos;x(:,1)]; ypos = [ypos;x(:,2)]; theta = [theta;x(:,3)]; v = [v;x(:,4)]; x_distance = xpos(end); if (x_distance > x_displacement_required) && (x_distance < x_displacement_required+x_displacement_accuraccy) achieved_distance = 2; current_monitor = [iteration omega x_distance]; monitor = [monitor;current_monitor];

elseif ((x_distance > x_displacement_required) && (omega_increment < minimum_omega_increment)) disp('Solution is not converging to the accuracy required.'); disp('Shown below are controlled parameters that may still be of interest'); achieved_distance = 2; current_monitor = [iteration omega x_distance]; monitor = [monitor;current_monitor]; elseif (x_distance > x_displacement_required) omega = omega - omega_increment; % This take omega back to a value that is definitely less than the required x displacement omega_increment = omega_increment/2; % Decrease the size of omega_increment so a more accurate answer is achieved end end solve_time = toc; tspan1 = tspan1.'; tspan2 = tspan2.'; tspan = [tspan1;tspan2(:,1)]; %% Calculate angles over time span phi_f_1 = a1*sin(omega*tspan1+Bf); % Front Platorm Angle phi_b_1 = a2*sin(omega*tspan1+Bb); % Back Platform Angle delta_1 = ar*sin(omega*tspan1+Br); % Rotor Angle

%% Calculate angles over time span phi_f_2 = a1*sin(omega*tspan2+Bf); % Front Platorm Angle phi_b_2 = a2*sin(omega*tspan2+Bb); % Back Platform Angle delta_2 = -ar*sin(omega*tspan2+Br); % Rotor Angle phi_f = [phi_f_1;phi_f_2(:,1)]; phi_b = [phi_b_1;phi_b_2(:,1)]; delta = [delta_1;delta_2(:,1)]; snakeboard = zeros(numel(xpos),8); snakeboard(:,1) snakeboard(:,2) snakeboard(:,3) snakeboard(:,4) snakeboard(:,5) snakeboard(:,6) snakeboard(:,7) snakeboard(:,8)

= = = = = = = =

tspan; xpos; ypos; theta; v; phi_f; phi_b; delta;

line_0 = ['User specified ',num2str(x_displacement_required),'m in ',num2str(duration),' seconds with an accuracy of ',num2str(x_displacement_accuraccy),'m.']; line_1 = ['For 0 to ',num2str(duration/2), ' seconds:']; line_2 = ['Then, from ',num2str(duration/2), ' to ',num2str(duration),' seconds:']; line_3 = ['Distance travelled: ',num2str(x_distance), ' (m)']; line_4 = ['Calculation time: ',num2str(solve_time), ' (s)']; a = '%4.12f'; disp(line_0); disp(line_1); fprintf('\t omega_f = '); disp(num2str(omega,a)); fprintf('\t omega_b = '); disp(num2str(omega,a));

fprintf('\t omega_r = '); disp(num2str(omega,a)); fprintf('\t a_f = '); disp(num2str(a1)); fprintf('\t a_b = '); disp(num2str(a2)); fprintf('\t a_r = '); disp(num2str(ar)); fprintf('\t B_f = '); disp(num2str(Bf)); fprintf('\t B_b = '); disp(num2str(Bb)); fprintf('\t B_r = '); disp(num2str(Br)); disp(line_2); fprintf('\t omega_f = '); disp(num2str(omega,a)); fprintf('\t omega_b = '); disp(num2str(omega,a)); fprintf('\t omega_r = '); disp(num2str(omega,a)); fprintf('\t a_f = '); disp(num2str(a1)); fprintf('\t a_b = '); disp(num2str(a2)); fprintf('\t a_r = -'); disp(num2str(ar)); fprintf('\t B_f = '); disp(num2str(Bf)); fprintf('\t B_b = '); disp(num2str(Bb)); fprintf('\t B_r = '); disp(num2str(Br)); disp(line_3); disp(line_4); plot(monitor(:,1),monitor(:,3),'s-'),xlabel('Iteration'),ylabel('x-displacement (m)');

13.

Stable Resonance Gait

% Title: Stable Resonance Gait % Filename:stable.m % Author: Jonathan Jamieson % Date: September 2011 – April 2012 % Purpose: Caclulates a stable resonance gait %% Clear Previous figures, variables etc clc; clear all; close all; min_a_p = 0.0; max_a_p = 0.8; min_a_r = 0.0; max_a_r = 0.7; duration = 1*2*pi; omega = 1; a_increment = 0.01;

time_increment = (2*pi)/omega; a_p_range = min_a_p:a_increment/10:max_a_p; a_r_range = min_a_r:a_increment*10:max_a_r; m = 75; % kg, mass of snakeboard and rider l = 0.285; % m, half the length of the snakeboard J = 0.22; % kg.m^2, moment of inertia of crossbar and platforms Jr = 14; % moment of inertia of rotor Jp = 0.013; % moment of inertia of platforms Bf = 0; Bb = 0; Br = 0; x0=[0;0;0;0];

time_step = 0.01; % The resolution of the differentiation current_time = 0; snakeboard = [0 0 0 0 0 0 0 0 0 0];

while current_time