Lyapunov Functions and stability

The intuitive definition of a stable system, is one where the system tends towards a stable equilibrium point. For a mechanical system this is best represented by a system which dissipates energy and comes to rest at an equilibrium point.

This is the basis of Lyapunov direct method, where the stability of the equilibria of a system are analysed using a Lyapunov function.

for a dynamical system:

we assume an equilibrium point at x = 0;

In order to determine the stability of the equilibrium point in terms of Lyapunov, we identify a differentiable function denoted V, as a candidate for the Lyapunov function. Fundamentally, we will determine from this function whether the system will remain at the stability point if initialised a it and whether the system tends towards said equilibrium point globally or locally.

for a subset defined on some region containing the region such that:

• V(0) = 0;
• V(x) > 0 for all with ;
• for all $$x \in S$

If the above statements are true that then x = 0 is a stable equilibrium point.

it is important to note that the Lyapunov function is not unique for a system and that is can be difficult or in some case in possible to identify a suitable Lyapunov candidate.

Where we have an autonomous system, such that that the dynamics of the system are not explicitly dependent on time, the Lyapunov function will not depend on tiem either. Therefore, the derivative of sai function can be defined as follows:

Local Stability:

function Again taking a function for which f(0) = 0 without loss of generality:

The system is locally stable about the equilibrium point where:

• V(x) is positive definite;
• is negative semi-definite.

if more precisely is negative definite, the equilibrium point is asymptotically stable.

Lyapunov theorem of global stability:

By definition a globally asymptotically stable equilibrium point is one that is stable, as defined above, everywhere i.e. for all . Thus:

• V(x) is positive definite everywhere in the state space;
• is negative semi-definite.

If the equilibrium point is globally asymptoticaly stable, it implies that it is the only stable equilibrium point in the system.

A third condition, allows us to identify if nto only the equilibrium point but the system itsels is globally asymptotically stable:

as i.e.

This means that V(x) tends to infinity as x tends towards positive infinity, in either the negative or positive direction. This has the effect of essentially bounding the bounding V(x).

Invarient Set:- A subset of the space $ such that if the system is initialised within the subset, it stays there and cannot leave. Once it is initialized in the subset it cannot leave the subset which comprises a region .

Example of an invariant set:

• The equilibrium point;
• The state space within the invariant set;
• any autonomous system’s path;
• limit cycle: a subset within the state space in which if the system enters it will follow a path in the state space.

Reasons for the conception of invariant set theorem:

• To provide a looser definition of asymptotic stability that is not constrained to a only systems with a negative definite Lyapunov function;
• Allowing the study of system convergence to more general behaviours, not just the equilibrium point.

conclusions are less precise and the hypothesis of the invariant set theorem are more general.

the Principle of the Theorem Let S, a subset in , be a bounded invariant set. Assuming there exists a differentiable function such that:

If we take B to be the largest invariant set contained within () the set of x values within the subset S for which the derivative of V(x) is equal to zero, then all trajectories starting in S approach the invariant set B as time tends to infinity. This means that the

According to LaSalle’s Theorem, for a dynamical system:

we no longer have to assume that Lyapunov function V is positive definite.

For a scalar function V we assume:

• everywhere in state space

• 

all trajectories tend toward the subset

More precisely, the system does not only tends toward R, but the largest invariant set within R.

Example:

where

0 \ for \ x \neq 0 “>

0 \ for \ \dot x \neq 0 “>

Therefore, the trajectory of both functions b(x) and c(x) have to pass through the point zero, as they both have the same sign as their argument.

In addition, the point and is a globally asymptotically stable equilibrium point.

virtual physics:

Assuming we are working on a mechanical engineering problem we can take the total energy of the system as the obvious Lyapunov function to consider.

We can see that the Lyapunov candidate above is scalar, however, the last term :

does not guarantee the following assumption

therefore we need to prove that we still need to show that

by taking the differential of the Lyapunov candidate (Potential energy equation)

This shows that , which corresponds to the set R.

For it to be globally asymptotically stable we said the system should tend toward he largest invariant set within R.

To verify that the system stays in R{} if it starts there, we can do this by computing the derivative which would need to be 0 for .

We can see from the initial dynamics of the equation that for then the function is also equal to zero, hence:

The only way the system tends towards an invariant set within R, is if c(x) is equal to zero, which is at x = 0.

Hence, the system is globally asymptotically convergent to the equilibrium point:

Important Note: This is a very simplified summary of my own notes which I personally find helpful as a shortlist of fundamental definitions. However, these concepts can seem rather abstract and complicated, and they can be at times, therefore to understand them better and really get a better grasp of Nonlinear dynamical systems I would advise reading through notes and books which cover these topics in more depth and provide practical examples which are extremely helpful in building a better understanding.

Velocity Kinematics

The velocity kinematics express the relationship between the velocity of the individual joints of the robot arm and the velocity of the end-effector in cartesian space.

The translational components of the end-effector pose, i.e. x and y component in 2D, for a 2 joint planar robot is expressed in previous posts in the form:

We can then simply differentiate these terms to fin the velocity components of the end_effector in x and y giving:

Rearranging the expressions of and above in terms of the joint angle velocities and , the expressions can be written in matrix as follows:

From this we can derive the relationship between the x and y components of the end-effector’s velocity an the velocity of the joint angles in the form of a Jacobian Matrix.

where $$\nu$$ is the end-effector velocity and the Jacobian is a function of the joint angles of the robot arm. This succinct form also applies to other robot arms with more joints and even in 3D space and can therefore be used to represent the relationship between the joint angles and end-effector velocity.

The Jacobian is the matrix equivalent of the derivative of a function f(x), where y = f(x). It’s the derivative of the function, who’s argument and result are a vector and the who’s elements are the partial derivative of each of the results with respect to the inputs to the function.

Based on the definition of of the end-effector velocity above, we can obtain an equation for the joint velocities required to achieve a set end-effector velocity.

In the 2D 2-link manipulator case the Jacobian J(q) was a 2×2 matrix, as there were only two joints to deal with.

However, there is one issue which may arise when trying to compute the velocity of the joint angles using the the inverse Jacobian, when the inverse Jacobian becomes essentially infinite due to a singularity in the system. This arises when joints align, such as in the above case, whereby the to joints essentially become one and the robot loses a degree of freedom. This is known as a singularity pose.

Mathematically, this is a result of the determinant being zero or close to zero, at a singularity pose.

On MATLAB we can check for a singularity in the inverse of the Jacobian by looking at the condition number :

>> cond(J(q));

Consider the set of all joint velocities with a unit norm.

Consider the end-effector velocity as a function of this set, such that the set of joint angles have a magnitude of 1.

Based on this we surmise that the set of velocities is a vector x vector product

This constraint is carried forward into the equation for calculating the required joint angle velocities such that:

 

This equation resembles an ellipse in the velocity space.

The ellipse, shown in the Figure above, is a function of the angle configuration, where the shape of the ellipse will vary depending on how the joint angles re configured.

The Ellipse demonstrates that, as per the position of the end-effector, the configuration will alter the manipulator’s ability to move it’s end-effector in a given direction and the speed at which it moves in said direction.

If we take a practical scenario where we intend to control a robot arm to move the end-effector along a straight line in cartesian space, we need to specify the cartesian velocity that the robot should follow and compute the required velocity of the individual joint angles.

However, a practical obstacle to overcome is the continuous change in joint angle configuration that results from this motion. The Jacobian which defines the relationship between the cartesian velocity of the end-effector and velocity of the joint angles is a function of the joint angles. Therefore, as the joints move the joint angles and therefore the Jacobian will no longer be correct.

Therefore, we need to periodically update the Jacobian by discretely updating joint angles and computing the Jacobian matrix at every time step, such that for an initial joint velocity:

we can compute a new set of joint velocities after every time step , where the new joint configuration is defined:

This process is repeated until the end-effector reaches its destination.

Velocity Kinematics in 3D

We use velocity kinematics to compute the the rate of change of the end-effector pose.

To describe the velocity kinematics o the end-effector we can use a numerical approximation akin to that of a finite element difference.

we can use this approach to compute the rate of change of the end-effector cartesian coordinates x, y and z as well as the rate of change of it’s orientation.

Over small changes, the approximation of a derivative can also be expressed
by taking a finite difference:

In the case of a robot manipulator we can apply this to the forward kinematics.

Where must be a very small change in the joint angle.

We can apply this to the end-effector pose to compute the approxiamte partial derivative which describes the rate of change of the pose (derivative) in terms of the initial and final end-effector pose :

The transformed pose can be written as a matrix, in terms of the rotational and translational components.

The derivative of the translational component of the transformation can then be used to compute the velocity of the end-effector in cartesian space, when multiplied by the change in the joint angle over time :

The velocity of the end-effector is defined:

For the Rotational component of the Matrix:

Where

S is the skew symmetric matrix and function of the angular velocity denoted

The skew symmetric or anti-symmetric matrix has the following properties:

• The Transpose is equal to the negative of the matrix:

• Is always singular, i.e. the determinant is equal 0;

• Any matrix is equal to the sum of a symmetric matrix and skew symmetric matrix;

In 3D space for a position vector of cartesian coordinates x, y and z. The skew symmetric matrix takes the form:



We can the vector cross products of a matrix vector product a x b = S(a)b

The derivative of the Rotational component of the Transformation matrix over time is defined as



where the , the above equation can be rearranged with the skew symmetric matrix as the subject:



where the Skew symmetric matrix takes the form:



From this we can show that for joint 1, the component of the angular velocity of the end-effector due to the rate of change of the first joint angle of the robot arm can be defined:



where a, b and c are arbitrary elements of the vector.

We can then describe the component of the transnational and angular velocity with respect to the rate of change of joint 1, as a column vector of six elements:



For a 6 joint robot such as the puma 560, we would then have a column for each of the joints meaning that the vector defining the relationship between the joint velocities and end-effector velocities would be a 6×6 matrix. As with the 2D case, this matrix which maps the joint velocities to the end-effector velocity is the Jacobian, which is a function of the joint angle configuration.



In MATLAB, using the Robotics Toolbox, we can compute the Jacobian for a vector of joint angles q simply using the method “jacobn”, e.g. p560.jacobn(q).

Provided it is not singular, we can invert the Jacobian to give us the expression for the joints velocities needed to achieve a specified spatial velocity at the end-effector.

Inverse kinematics

Inverse kinematics is the function that gives us the joint angles needed to achieve a particular end-effector pose.

This is described:

To explain how we can do this, we will look at the inverse kinematics of a simple 2-joint robot arm for starters.

The objective is to calculate the joint angles and needed to achieve the end-effector pose

It should be noted that either or can be defined in terms of the other.

We start by defining r, simply using Pythagoras i.e.

we then need to use the cosine rule to determine the angle alpha, between the two links.

cosine rule:

Where, b and c are represented by the link lengths and , respectively, and length a is r calculated above. From this we can compute an expression for angle :

which we can see from Figure 1 is also related to angle by . Therefore, considering that , we can compute the joint angle to be:

Taking a second triangle about the second link and joint angle , we can compute the sides of the triangle in terms of the joint angle and link length :

We can then use these lengths to find the angle with the relationship:

Given the sum of and joint angle relate to the triangle of base x and height y, shown below, we can see that the angle Beta relates to the joint angle by:

We can then substitute the expression for to compute the equation of in terms of the link lengths and second joint angle :

Looking at our expression for we can see that our inverse kinematics give us two possible solutions for the joint angle, positive or negative.

Either one will give us the desired end-effector pose, as the expression for joint angle is itself a function of . However, they correspond to different joint angle configurations.

The method covered above is a geometry based approach. However, we can solve the inverse kinematics taking a more analytic approach, using algebra.

For the same 2-joint robot arm we discussed in the previous example, shown in Figure 1, we have an end-effector pose that can be described by the forward kinematics expression:

This is a sequence of two translations and two rotations.

Multiplying out and expanding this expression, gives a homogeneous transformation matrix of the form:

In this case, what we are most interested in is the position of the end-effector defined by the x and y coordinates, which are described by the 3rd column of the matrix above.

We can add the square of the expressions defining the x and y coordinate of the robot end-effector, which resolves to:

we can rearrange the equation to isolate the term and give us an expression for the required joint angle in terms of the end-effector coordinates:

This is the same solution computed using the geometry based approach.

Now to find the solution to . We can take simplify the x and y coordinate expressions substituting short-hand notations:

This gives us:

These expressions are of the form:

which has the solution:

Therefore, by substituting the variables of the expression y into the identity above such that , and we can solve for :

We can see this mirrors the solution obtained from the geometry based approach.

The inverse tan function:

sometimes also denoted atan or arctan suffers from a limitation.

The standard inverse tan function is incapable of recognizing an angle passing into the (LHS) or imaginary plane. this is because the standard tan function only has range and therefore cannot account for a rotation in the LHS or a negative sign. The inverse tan function atan is used to get around this issue, which has a range and accounts for the sign of the coordinates input.

The function atan2 is found in many mathematical libraries, including MATLAB, where it considers the sign of both coordinates to evaluate the angle:

In MATLAB, for a modular robot R, we can use the command ikine in the Robotics Toolbox to compute the joint angles required to achieve the specified end-effector pose described by a transform matrix denoted T, as follows:

The R is an arbitrary notation I have assigned for a serial link robot, but this needs to be defined in the MATLAB library. To simply test the command for the time being, you could use the mdl_hyper3d() or mdl_puma560 which are already have defined libraries in the Robotics Toolbox. In which case the for the hyper3d, the inverse kinematics can be computed using the command:

Where there are more than one set of joint configurations which achieve the desired end-effector pose, as shown above with the positive or negative values of , the ikine command allows us to specify the desired configuration: ‘r’ for right and ‘l’ for left.

We have covered the analytical method of computing the joint angles using inverse kinematics. However, it can be more convenient and practical to use what is called the Numerical method, especially when the robot has a large number of joints, i.e. above 6.

The advantage of using the analytical approach is that it explicitly shows the multiple configurations that the robot arm can be in for a particular end-effector pose. When working with a robot with a small number of joints, the analytical approach can be the best option as it provides a compact solution which is faster to execute than the numerical alternative.

However, it may be more mathematically difficult or even impossible to find the inverse kinematics of a more complex robot. In this case the Numerical method should be used.

We have seen that we can always find the forward kinematics of the robot. To achieve a desired pose we need only adjust the joint angles so as to move our end-effector from the starting pose to the final pose we want. This can be presented as an optimization problem defined:

Using MATLAB, we use the ikine command mentioned above. This is fairly simple, and requires us to only specify the final end-effector pose and the initial joint angles, as shown by the example below for a puma560 robot:

One shortcoming of the Numerical method, is that it can be difficult to obtain all the possible joint angle configurations which can achieve a particular end-effector pose, as MATLAB will only return the one possible configuration based on the initial joint angles specified.

In addition, if the initial joint coordinates are too far from the joint coordinates needed to achieve the desired pose, than the solution may not converge.

Note: When trying to compute the inverse kinematics of a redundant robot we need to use the numerical solution. But there will be a near infinite number of joint configurations that result in the same end-effector pose.

Singularity is an problem in robot control similar to Gimbal Lock, which rises when two or more joint axis are aligned.

The alignment of the two joints results in a loss of a degree of freedom as the motion of either joint simply creates the same motion in the end-effector. It, therefore, becomes impossible for the robot end-effector to achieve certain orientations while these joints are aligned, that were otherwise possible when the joints where not aligned.

Due to this singularity, the robot will be unable to adopt certain orientations that were otherwise possible when the joints where not aligned. We will see this problem pop up in several robotics problems.

Moving a robot from one pose to another, we need to determine the joint angles at each pose, by carrying out the inverse kinematics of the poses.

To move smoothly between the two poses, or set of joint angles, we need to interpolate the vector quantities overtime.

We can simulate the motion between two end-effector poses, within the robot workspace, by creating a trajectory between the set of joint angles, and specifying the number of time steps in which to perform the motion. We can use the command jtraj:

In this example I have asked MATLAB to perform the notion over 50 time steps, between a set of joint angles q0 and q1, which will return a matrix of 50 rows and 6 columns (for 6 joint robot in 3D such as the puma560), where each row corresponds to a time step and the columns correspond to a joint angle.

We can then visualize this using the plot command with the robot we are considering i.e. for a pum560:

We can then use the qplot(Tr) command to plot a figure of the change in joint angles over time. This should show that the Cartesian position of the end-effector moves smoothly as a function of time, as we would want.

However, if we plot the poll-pitch-yaw angles against the time steps. We can see from this that not all the angles will remain at zero constant. This is because we have only constrained the orientation at the beginning and end of the trajectory, at poses A and B.

This indicates that the end-effector did not follow the desired straight path and that its orientation has deviated along the trajectory.

We can use Cartesian Interpolation to resolve this.

As opposed to interpolating the joint angles, in cartesian interpolation we will directly interpolate the Cartesian pose and the perform the inverse kinematics to find the joint angles at every individual time step between the two specified poses.

where s varies smoothly from 0 to 1.

So in Matlab we will use the command ctraj for the cartesian trajectory from an initial and final pose and pass in the number of time steps, e.g.:

This will give a series of 50 poses, each of which are 4×4 matrices, one matrix for every time step.

We can the pass the cartesian trajectory defined (Tc) into the inverse kinematics of the robot (puma560):

This will return a 50 sets of joint angles, one set for every time step.

Configuration space – the space in which the joint coordinates of the robot are described. Robot configuration is described by a vector of generalised joint coordinates, denoted:

Where N is the number of joints in the robot and configuration space C is the set of all possible joint configurations.

The generalised joint coordinates can be either angles (revolute joint) or lengths (prismatic joint).

The dimension of the configuration space is given by the number of joints in the robot.

The task space defines the space of all possible end-effector poses. The task space is the space for all possible end-effector poses and the task space number defines the degrees of freedom of the robot.

The dimension of the task space is the number of degrees of freedom of the robot or the number of elements which define the pose of the end-effector. This can be a maximum 6 (degrees of freedom), .

In robotics, a robot with a configuration space dimension far larger than the task space is known as a redundant robot. This is a robot capable of achieving any end-effector position and orientation within its workspace.

kinematics is a branch of mathematics that studies the motion (position and velocity) of a body, or system of bodies.

In Robotics there are two types of kinematics problems to consider:

-Forward kinematics;

-Inverse kinematics.

The forward kinematics of a system, e.g. a robot arm, give us the pose and coordinate of the robot end-effector given a set of joint angles.

Consider a simple single joint robot arm.

Take E to be the pose of the end-effector which is has a translational and rotation component such that:

In the 2D case, for an angle and link length , the pose of the end effector would be expressed:

From this we can see that the coordinates of the end-effector are represented by the last column of the matrix:

The command in the Matlab robotics library to model the forward kinematics of the single planar system is:

where q is the angle of rotation, in this case q1.

Tip: To visualize a representation of a robot arm, we can import the model of a planar 1-joint robot using the command from the Robotics Toolbox, and input the joint angles using the command ‘‘.

This way of representing the forward kinematics of the robot arm translates identically to a 2 or 3 link robot arm, where we can express the end-effector pose as the sum of the rotations and translations of each link of the robot:

In Matlab this would be written using the ‘trchain2()’ command as follows:

Once again the translational component of the homogeneous matrix represents the x and y coordinate of the end-effector which is defined by the third column of the 3×3 matrix of the pose. In the case of the simple two link robot arm:

We know that for a two planar robot such as the one shown in the figure above there are always two joint angle configuration that will position the end -effector at any given coordinate within the workspace.

With a 3 joint robot arm, the end effector of the robot can reach any position within its workspace and achieve any orientation for that given position.

Take the simple 3-link robot arm shown below.

we can describe the pose of the robots arms end-effector as a chain of homogeneous transformation matrices which belong to the Euclidean special group of 3 dimensions.

A very commonly used industrial robot which which uses both revolute and prismatic joints is the Scara robot, shown below.

and translation. For example, when computing the motion of the prismatic joint in Matlab, using the trchain command, it is important that the correct axis of translation is defined. In the previous examples we only covered 2-DOF arms which had only revolute joints rotating on the z axis, therefore, we only input translations in the x axis and did not need to specify the axis of the rotation matrices. However, in the case of the Scar robot arm, the end effector is attached to a prismatic joint and the robot end effector can move in 3 dimensions. As such the End effector pose would be expressed:

The Denavit and Hartenberg system was developed as a general theory to describe and articulated sequence of joints, where each joint is defined by four parameters.

link 0 is used to describe the base of the robot itself. The following links will then be numbered 1, 2…n from the base to the end-effector.

The principle of the DH (Denavit and Hartenberg) method is that a coordinate frame is assigned to the end of each link (at the joint or end-effector) and describe the pose of said link frame with respect to that of the previous link frame.

As usual, the relative pose from one link frame to the next is represented by a homogeneous transformation matrix of the form:

It is important when modelling robots arms with 6-DOF and above to specify the axis of rotation

As we can see from the expression above, the transform from link frame i-1 to link frame i comprises two rotations and two translations: a rotation and translation about the z-axis, followed by a translation then rotation about the x-axis.

From this we can describe the relationship between the two frames in terms of just 4 parameters, in this case ( and )

Pose has a translational component and rotational component (euler angles or roll-pitch-yaw). Therefore, a transformation in 3D would have 6 parameters.

Using the Denavit Hartenberg notation, imposes certain constraints on the position of the link frames:

• The x axis of frame i intersects the z-axis of i-1;

• The x-axis of frame i is perpendicular to the z axis of i-1.

These constraints allow us to describe the change in relative pose using only 4 parameters.

For a robot arm we can simply stack groups of these link transforms (each defined by 4 parameters) to compute the end-effector pose relative to the robot base, link 0.

Normally only one or two terms of the link transforms are variable or changing, and the other parameters are typically constants as they are simply a function of the mechanical design of the robot.

Take, for example a 3 link manipulator shown in Figure 4. In this case, only the rotation about the z-axis is changing, while the the remain three terms remain constant for each link.

This will obviously vary based on the design of the robot, for example for a prismatic joint the translation along the z-axis becomes the variable parameter.

The Denavit and Hartenberg notation can quite easily be computed in Matlab, using the Robotics toolbox.

By declaring a a simple 4xn matrix to define the four parameters of the N number of links. A 4×2 matrix will define a 2 link robot arm, where each row represents a link.

We can then use the Toolbox function SerialLink(DH) to create a serial link object of a robot arm. An example of this is displayed in the Figure below.

There are two ways to model the system in Simulink. Using the standard Simulink library and the Simscape library. They differ in several respects as you will see in the following section.

The challenge in modelling an inverted pendulum in Simulink, is the physical constraint created by the pin joint between the cart and pendulum, which effectively reduces the degrees of freedom of the system to 1DOF each, the position of the cart and , respectively.

Once again taking an Inverted Pendulum system as illustrated in the Figure below.

We can use Newton’s Second law to generate the differential equations for those degrees of freedom such that:

Where b is the friction coefficient and N the reaction force.

The sum of the moments acting on the Pendulum, in terms of the vertical and horizontal reaction forces P and N, is given:

These two equations can be rearranged to give expressions for the acceleration of the cart, resulting from this external impulsive force F, and resultant angular acceleration of the Pendulum.

In order to model the dynamics of the system in Simulink, we need to consider the interactions of the forces between the cart and pendulum P and N by defining the two reaction forces.

We previously need to mathematically derive the equation of the system and resolve the resulting algebra in order to model the system in Matlab. However, when using Simulink, we can rely on Simulink to resolve the algebra for us, asking this approach less mathematically taxing.

As the coordinates of and are exact functions of we can represent their derivatives in terms of .

For the x-component:

The y-component:

The equations for the derivatives of the x and y components of the pendulum position can the be substituted into the equations of the reaction forces P and N to give the following:

We can use the reaction forces and accelerations above to define the dynamics of the system.

As you can see this approach, using Simulink, can be far less algebraically taxing than computing a transfer function or state-space form for the system and can be very convenient when testing a model.

We the need to build our Simulink model blocks. Starting with the derived differential equations for for acceleration of the cart and angular acceleration of the pendulum and the two reaction forces, we need to insert four function blocks one for each: , P and N.

When defining the equation of each of the functions, we need to replace the variables (no-constants) as inputs e.g. u(1), u(2) etc.

For example can be written in the form: (1/M)* (u(1)-u(2)-b*u(3))

where inputs u(1) is the impulsive force F, u(2) is the reaction force N and u(3) is the velocity , which is the integral of this same block in the previous time step.

Then we can use the integrator block to compute the and x from and and from the function block.

The we can proceed to pass the and functions through their two integral blocks and connect the the variables to their desired input ports as shown indicatively n the figure below.

It is important to be careful when connecting the blocks, especially where the function blocks have multiple inputs, as the order in which the inputs are connected can fundamentally change the expression of the function. The number of the input port corresponds to the input number specified in the function’s expression i.e. u(1) represents the first input port.

We can the add an input block F ad two output blocks, x and Theta.

Figure below shows an indicative layout of what this might look like.

Note: The layout of the links will vary depending on how you have written the function expressions.

Once these blocks are in place, and inputs correctly linked to their respective sources, we can turn this into a subsystem, by simply highlighting all the blocks, right clicking and selecting ‘create subsystem from selection’ from the drop-down menu.

We can then test the response of the subsystem to an impulsive force on the cart. One way to achieve this is by connecting an pulse generator block to the subsystem input to simulate the impulsive force and a scope block to the output ports to measure the cart position and pendulum rotation, as shown in the figure below.

When using the Simulink library we will still need to input the differential equations for the acceleration and reaction forces derived above. However, the Simscape method is a more literal translation of the system, making use of joints and virtual actuators. The Simscape approach has the advantage of giving us a visualization tool which can provide a visual representation of the system response, once we have run the simulation.

In this section we will cover how we can use the Simscape library extension of Simulink to model the nonlinear inverted Pendulum system. This involves the least amount of algebra of all the approaches as blocks in the Simscape library represent actual physical components. As such, we can build complex multi-body dynamic models without the need to build mathematical equations from basic principles as per the original Simulink case. In addition, this method benefits from the fact that using the blocks from the Simscape library allows us to produce a visual representation of the response of the system once we run the system.

Building the system using the Simscape library will require:

• Two body blocks
• 3 joint actuators;
• Two Prismatic joint block;
• One Revolute block;
• Two Ground Blocks;
• Two joint Sensors;
• A Gain block;
• One Machine environment block;
• Two outputs and an input;
• One Joint Initial Condition block.

Within the Body block parameters we must also define the location of the centre of mass and the shape and dimension of the bodies by specifying the coordinates of each corner of the body, with respect to the centre of mass specified. For a 2D case, we only need to define 4 corners of the cart, the coordinate of the centre of mass, and specify what the coordinate frame used, either world or CG (centre of gravity) to define the position of these edges.

In addition to defining the body parameters, two further inputs are required to define the position at which the impulsive force and frictional force act on the cart. An example of the body block parameters for the cart and pendulum are shown below.

For the ‘Port Side’, we simply assign the centre of mass and the two force inputs, to the ‘Left side’ and the remaining entries to the ‘Right side’.

<p

We can simulate an impulsive force using the pulse generator block. Once the amplitude and period of the pulse is set, the block should be connected to the input of our subsystem.

As per the Matlab approach, we can linearise the system created in Simulink allowing Simulink to resolve the algebra for us. This steps to linearise the system are the same whether using the simscape library or standard Simulink libraries.

Step 1: Set the input signal ‘Force’ to the subsystem to ‘Open-loop input’ in the linear Analysis Points on the drop down menu; Set the two output signals ‘Theta’ and ‘x’ to ‘Open-loop Output, as per above.

Step 2: Select Linear analysis in the Control Design drop down menu under the Analysis toolbar. This will automatically prompt open the linear analysis window, as shown below.

To start with, select the Trim model function under the ‘Operating Point’ menu to create a new operating point. As per default, we should find that only the ‘steady state’ state specifications are ticked, and the maximum and minimum values for each of the states are set to infinity and negative infinity. We can leave this a it is and Start Trimming. This will create the new operating point.

Step 4: Select the created operating point in the Linear analysis workspace, then run the Impulse Plot in the Linear toolbar. This should automatically linearise the response of the system and produce two graphs, one for the angular deviation (theta) of the Pendulum and the other for the displacement of the cart, resulting from the impulsive force on the cart simulated.

Modelling an Inverted Pendulum

There are a number of ways to model and control an inverted Pendulum using MATLAB and Simulink. Each approach has its benefits. Four approaches will be outlined in the following posts, two in MATLAB and two in Simulink.

The modelling of an inverted Pendulum is a useful exercise to familiarise oneself with MATLAB and Simulink while applying your knowledge of dynamics. It is a particularly good exercise since the inverted pendulum problem is a practical practical problem which closely relates to a number of common place engineering problems.

We will model the inverted pendulum system directly using the MATLAB command window, by first deriving the necessary equations with which we can define the system and then inserting the equation/s into MATLAB, to allow the software to model the system for us. This is obviously far more math heavy than the Simulink approaches we will cover later on.

Two ways of representing the motion of the pendulum and cart in MATLAB are as transfer functions or in a state-space form. In either case we will need to look at the dynamics of the system and resolve the forces acting on the Pendulum and cart, for an impulsive external force F. A simple inverted Pendulum attached to a wheeled cart is shown below:

From the Figure above we see that the sum of the horizontal forces acting on the cart give the following equation of motion:

Where b is the coefficient of friction on the between the car and the ground, and N is the reaction force on the cart. For now we will only look at the sum of the horizontal forces.

Similarly, we then sum up the forces acting on the pendulum to show that the reaction force N for the pendulum can be defined by the angular velocity the equation:

By substituting the equation for N back into the original equation of motion of the cart, we obtain:

The sum of the forces perpendicular to the pendulum can be described by

In order to substitute out the reaction force terms P and N, we can simply use the equation for the sum of moments about the Pendulum, which is a function of the reaction forces and the moment of inertia I:

By substituting this into the equation describing the forces about the pendulum, above, we can rearrange to get:

As we are, at least for the purpose of this exercise, only applying this control design to a linear system, we can linearize the equations, such that , where the defines the deviation of the pendulum position from the upright equilibrium position. This is a reasonably valid assumption provided that under the control law we implement, the pendulum does not deviate any more than 20 degrees.

On this basis, we can also then linearize the trigonometric terms in the equation such that:

The linearized forms of the equations of motion derived will give the following:

and

Writing the Transfer Function of the pendulum and cart:

In order to compute the transfer function describing the linearized system, we must take the Laplace transform of the equations, assuming zero initial conditions, such that u(0)=0, x(0)=0 and .

The resulting Laplace transforms of the two equations are shown below:

and

Note: A transform function represents the relationship between a single input and single output at a time. Therefore, we will need to compute a separate transfer function for the position of the cart and the angular deviation of the pendulum.

We will then need to rearrange the derived equations in terms of it’s inputs and outputs such that the Transfer function is defined:

The transfer function of the pendulum will be defined by the input and output . Therefore, we will need to substitute out the term of the cart position in order to compute the transfer function in terms of a single input, rotation angle . We can use the previous two equations derived to obtain an expression of X(s) in terms of .

This expression can then be used to substitute out the term in Laplace transform of the equation of motion to obtain an equation containing only a one input term and output term :

We can then rearrange this equation to obtain an expression of the transfer function of the form:

The transfer function equation contains a large number of constants (M, m, L and I) which substantially lengthen the transfer function expression and increase the possibility for error. For our own ease of use, we can shorten the transfer function by simply substituting out some of these terms with a constant q.

For example we can take

By substituting out these terms for q, we get a less chaotic expression:

We can see that in this transfer function, for the motion of the Pendulum, there are common real zeros and poles. This can be simply canceled out to give:

This is the final transfer function for modelling Pendulum orientation. When writing this in MATLAB, it is important to remember to define the constant ‘q’, use to simplify our derivation, in the workspace along with the other constants of the system before running the simulation.

In addition to the Transfer function of the pendulum, we can proceed to simply formulate the transfer function of the cart, to give cart position as the output, using the above expression by simply substituting back in the variable in terms of X(s). this will give the following:

We can now start transferring our work to MATLAB. We will need to start by defining our known constants (m, M, L, g, I, b and q) in the workspace, before transferring across the two transfer function expressions.

The steps for modelling the system in MATLAB are outlined further in the ‘Modelling in MATLAB’ Section below.

The command for modelling a system based on a transfer functions is:

> sys \_ tf(Tpend, Tcart)”>

We can also represent the model in state space form. To do this we need to identify the states with which we can define the system and compute the representative matrix. This method may seem more direct at first as it involves fewer steps and derivations but can be prove to be more mathematically tedious than the transfer function method, as it may result in having to work out more complex matrices, increasing the possibility for error.

Definitions: State space for is the Euclidean space in which the variable on the axis are the state variables, the state of the system can be represented a a vector within the space. As a result, when expressing the dynamics of the system in state-space form, the system in the inputs, outputs and states are also expressed as vectors.

We assign a state to each of the elements of the system’s dynamics and look at their dynamics separately.

In this case, using the two equations of motion previously derived for the inverted pendulum, stated above, we can identify 4 principal states that we can use to define the behaviour of the system:

We first need to define these states using the two equations of the system derived from the sum of forces acting on the cart and pendulum. It is important to note that in this instance we are still assuming a linearised system whereby the following is still true:

We can use the two above equations to formulate expressions for the four states. The states and can simply remain as they are. However, the other two states will need to be defined in terms of the former two states and the input u.

It is important to note that each of the state expressions can only be defined in terms of constants and lower order derivatives i.e. cannot be defined in term of and vice versa. Using the two equations above we can substitute out the unwanted variables from the respective expressions, and with a bit of rearranging, show that:

From this we can write the state space equation in matrix form:

The output of the state-space can be expressed in the matrix form as follows:

Where the first and second rows of the output are the position of the cart ad angle of rotation of the pendulum.

To transfer this over to MATLAB, we first need to define all our known constants, variables and matrices. For simplicity, we can express the state space equation in the form:

and

So as not to overburden ourselves when writing in the expression.

We can then define our states, inputs an outputs. To model the system using the state-space form we can use the command:

>”>sys_ss=ss (A,B,C,D, ‘statename’, states, ‘inputname’, inputs. ‘outputname’, outputs)

This state space form of the system can also be converted into a transfer function using the command: >>sys_tf=tf(sys_ss)

In order to begin modelling the system in MATLAB, no matter which approach is taken, after the known constants and variables of the system have been defined in the workspace.

We then need to define our inputs and outputs as shown below:

>input=\{‘u’\}”>

>output=\{‘x’,’phi’\}”>

We can then use the relevant system command, of the the respective type of expression, to model our system i.e. sys_tf or sys_ss for the transfer function and state space-form expression, respectively.

We can set the names for the input and outputs using the commands:

>set(sys_tf, ‘inputname’, inputs)”>

>set(sys\_tf, ‘outputname’, outputs)”>

In the specific case of the state space method we will also need to define our 4 states before inserting the system command:

>states=\{‘x’, ‘xdot’, ‘phi’, ‘phidot’\}”>

The Matlab script used for to model either cases are presented below. However, this could always be written in a different format. I have also added the script you can use to produce graphs showing the model’s response to an arbitrary impulsive response.

Script for modelling the inverted Pendulum using transfer functions in the Laplace domain:

Script for modelling the inverted Pendulum in the state space-form is shown below:

We ca use the same impulsive function and the script as shown in the ‘transfer function’ code to test the function for an impulsive force and produce graphical representation of the response of the cart’s movement and angular deviation of the Pendulum.

• trajectory is Defined as a path and a schedule, where as a path is not dependent of on a schedule (time)

There are number of ways of computing or projecting the trajectory of a system these will be discussed and compared below.

Quintic polynomials are commonly used in robotic for projecting trajectory. A quintic polynomial is a polynomial of the 5th grade, comprising 6 coefficients:

In addition to six coefficients, the polynomial has six boundary coefficients at time: t=0 and t=T:

These define the initial an final position, velocity and acceleration of the system. The derivatives of the displacement which is a function of time are:

By substituting in the initial boundary condition at time t=0, we see that:

For time t=T, the 6 constrained boundary coefficients can be described by the following 6×6 matrix:

We can plot the polynomial trajectory in 1D using the matlab function ‘tpoly()’ from the robotics toolbox:

This polynomial trajectroy plot gives us a smooth continuous curve of displacement velocity and acceleration as shown indicatively by the plots below.

From the second graph, displaying the variation of we can see that the system’s average velocity over the set time step is restricted to no more than approximately 50% of the maximum motor speed, meaning that the object is travelling much slower than it theoretically could be, therefore, taking longer to reach it’s destination than needed if it were using it’s maximum speed. From a practical robotics standpoint, this results in a waste of motor performance.

The trapezoidal trajectory is an alternative to simulating trajectory using the polynomial trajectory method. In a trapezoidal trajectory the system accelerates to the peak velocity, where it remains at constant velocity for as long a necessary before decelerating as quick as possible in order to reach it’s destination. This results in a trapezoidal shaped velocity profile, from which the profile get’s its name.

The displacement in the trapezoidal profile can be defined by the expression:

The profile comprises three main phases, a shown in Figure 2 below:

• the coasting phase – where the system has a constant velocity, represented by a straight line in the profile;

the acceleration phase;

the deceleration phase;



The area under the velocity profile represents the distance traveled. While this is be the same for the polynomial and trapezoidal trajectories, the time units required to cover that distance is significantly less in the latter.

The Trapezoidal profile is also referred to as a linear segment with parabolic blend given the nature of the position and velocity profile in these phases. In the acceleration and deceleration phases, the velocity is a linear function of time, meaning position in this phase is a parabolic function off time. The position and velocity profile in the acceleration and deceleration phases are defined by expressions:





In the coasting phase, the velocity is constant i.e. . This phase is represented by a straight line segment in the velocity profile.

The trapezoidal trajectory is represented in Matlab by the function ‘lspb’ in the Robotics Toolbox, which allows you to plot a trapezoidal profile by specifying the starting position, final position and time step:

>lspb(q1, q2, time \ step)”>

The function also allows the user to specify coasting velocity of the system if desired. However, it should be noted that the range of velocity can be quite restricted by the minimum and maximum velocity required to achieve the desired displacement.

The main advantage of the trapezoidal profile over the polynomial trajectory is that the system accelerates to a coasting velocity as quick a possible, allowing it to reach it’s destination much quicker that it would have following a polynomial trajectory.

The main disadvantage of the trapezoidal profile is that acceleration is not continuous. This, creates discontinuity which can be a problem for some robots. However, we often depend on the low bandwidth of the motion control system to smoothen out this discontinuity.



We can also plot a trajectory by using specified Intermediate points, also called via points. The robot can then move from one point to the next, in segments.

However, in performing these motions, the velocity of the robot will not be smooth, as the velocity is increasing or decreasing abruptly about each of the specified via points. An example of a trajectory with via points is shown in the figure below.



We introduce what are known as blends to resolve this by introducing a smooth mathematical curve to smooth out the velocity discontinuity around the via points. Once we get near a via point we transition the velocity of an acceleration interval, over a set period.

The profile of the trajectory consists of segments of linear velocity between points, with higher order or parabolic polynomial blends joining them together.

We see that the trajectory of the system moves from one to the next between the linear segments. However, we should note that the robot is not actually touching the via points, but instead miss them, simply passing near them before accelerating towards the next via point. This is a trade off between accuracy and speed. Depending on the acceleration time specified, the velocity will need to change to transition between the points resulting in larger

To compute a trajectory with specified via points, we can use the function ‘mstraj’ from the robotics toolbox, where we can: specify a set of via points, set the velocity, specify the time of each segment, input the starting point, specify the time interval, set the acceleration time. However, one of the variables should be left unspecified i.e. velocity , time specification or acceleration, otherwise the function will have too many input arguments to compute a feasible trajectory.

Increasing the acceleration time will result in a more rounded trajectory, and the distance by which the trajectory misses the specified via points increases.



Linear interpolation is a useful tool to predicting trajectory in a smooth and in a coordinated fashion.

The fundamental formula for a the linear interpolation of a position as a function of s is given as follows:



From this we can see that:

• x(0) is the initial value of position, the initial boundary condition;
• x(1) is the final value of position;
• x(0.5) is a half-way point within the set of possible position bounded by ;
• x can be a vector i.e. of the form (x, y, z);
• s is a scalar which can be between the set 0 and 1.

given that x is a function of s, it can be a smooth function of time, a per the polynomial and trapezoidal trajectories.

In Matlab, we can use the function jtraj from the workspace to linearly interpolate a half way point between set values and where these values can either be defined a vectors or scalars. This can be used as:



The jtraj function also allows us to specify the initial and final velocity of the system. We should be careful to ensure that the number of expressions for the velocity term agree with that of the position i.e. 2 and 3 elements in the 2D and 3D cases, respectively.

It is important to note that while this formula, can be applied directly for a set of vectors or even matrices by substituting for the function x, this formula cannot be applied directly to rotation matrices. This is due to the orthogonal nature of the rotation function, where each column of the rotation matrix is a unit vector. As the sum of two unit vectors can never be a unit vector, the sum of two orthogonal matrices can never be an orthogonal matrix.

Instead, to represent rotation using linear interpolation, we can simply use vectors to represent the orientation by translating the rotation matrix to find an initial and final angle.

By introducing the initial and final vector angles, we are able to apply the linear interpolation formula directly as follows, where the denotes the vector of angles, to interpolate any intermediate angles and plot a trajectory:



Where andn are the initial an final angle of the rotation.

We can do this by converting the matrix which describes the rotation into an initial and final set of angles. We can represent the set of angles as a vector of Euler angles or Roll-Pitch-Yaw angles.

As each element (angle) of this vector can be considered a point on a circle, where the rotation from the initial angle to the final angle is represents a segment along the circles perimeter, we can see that there are always two paths from the initial to the final angle. Given the convention to measure angles in the counter clockwise direction, computing this rotation may result in us following the longer of the two paths represented below:

Figure 4: Rotation Path Diagram

It is important to be aware of this so as to ensure that we always take the shorter of the two paths.



Quaternion interpolation, also known as ‘Spherical Linear interpolation’ is the most direct method of plotting the trajectory between any two orientations.

The formula for interpolating between two quaternions, which represent the two given orientations, is given:



The angle is derived from the elements of the initial and final quaternion such that:



This interpolation corresponds to a constant angular velocity about a fixed axis in space.

Quaternion interpolation can be performed in Matlab using the function



Where q2 and q2 are defined unit quaternions in the Matlab works-pace.



Poses comprise a translational and rotational component. In order to interpolate a pose, the two components have to be interpolated separately.

consider two poses denoted and

We can convert the rotation matrix of each pose into a unit quaternion given as follows: and

and interpolate between the two quaternions using the quaternion interpolation formula above.

The transnational components of the poses can then be interpolated as vectors using the basic linear interpolation formula previously discussed:



As discussed, can simply convert the interpolated quaternions back to a rotation matrix and combine it with the interpolated translation to compute a homogeneous transformation function.

Note: The Matlab commands used to model the above transformations and produce the graphics above, were not taken from a personal library, but make use of the ‘Robotics Toolbox for MatLab version 10.x’ library provided by Peter Corke, a professor of Robotic vision at QUT (Queensland University of Technology) and Director of the Australian Centre for Robotic Vision (ACRV) . The Toolbox provides many functions useful for studying and simulating classical arm-type robotics, such as kinematics, dynamics, and trajectory generation. I have personally found the functions available in the robotics toolbox very useful in my studies and in modelling robotic systems in Matlab and I will regularly make use of it and reference the Toolbox throughout my future posts. A link to the Robotics Toolbox can be found on the QUT Robot Academy website, which also covers a lot of the fundamentals of robotics in a masterclass that I think may be helpful and practical introduction to Robotics which I highly recommend.

Link: https://robotacademy.net.au/resource/robotics-toolbox/

Robotics Toolbox

Quaternion Representation of Rotations in 3D

A Quaternion is a hyper complex number. An ordinary complex number would comprise a real and imaginary part of the form: ; where b represents the imaginary part of the expression and i is the imaginary term. The imaginary term of a complex number can often be expressed by it’s square .

A hyper complex number, is one that comprises one real part and three imaginary parts, and therefore has three imaginary terms. The hyper-complex number is the basis of the quaternion and the fundamental formula of a quaternion is described by the expression:

The number can also be written in the form of a scalar (real part) and a vector, which represents the three imaginary parts –

This is also commonly written in the form:

Quaternions can also be represented by 2×2 complex matrices.

A unit quaternion is a quaternion who’s magnitude is equal to +1,

Unit Quaternions are used to encode rotations in 3D space. They are similar to the angle-axis representation previously covered in that respect, however, using quaternions to predict the trajectory of a system and interpolating a point in a rotating or translating system is more accurate. This will be discussed further in the following section.

The real part of the quaternion is defined: ; The complex part is defined as: ; where denotes the rotation axis.

Quaternions are non commutative but can be compounded, as shown below:

The inverse of the quaternion can also be comuted usign the Hamiltonian rule:

In Matlab, we can simply use the function “UnitQuaternion(R)” from the Robotics Toolbox, to compute the Unit Quaternion of a rotation R, which we have defined in our workspace.

Coordinates

To Describe rotations in the 2D plane we can use 2 Coordinates frames A and B, shown below:

The Vector to Point P can be defined by the following expressions, in terms of either frame of reference:

Where, and represent the distance relative to the frame of reference A and B. The values and can be expressed in terms of the frame of reference {A}:

Substituting these two expressions into the expression of vector and rearranging the terms as coefficients of and gives the following:

Given vector and are equal, and can be expressed:

This can also be written in Matrix form:

This can also be written i short form, where R is the rotation matrix.

There are a couple of rules to be aware of when resolving rotation matrices:

• The inverse of a rotation matrix is equal to the transpose of the matrix e.i.
• The determinant of the rotation matrix is equal to positive 1, i.e.
• Rotation matrices in the 2D, belong to the Special Orthogonal group of dimension 2,

Rotation and Translation in 2D

Take the following scenario, where we have two frames of reference {A} and {B}, where the frame are separated by translation denoted by vector V and Rotation R. For a point P, the vector to point P in terms of the frame of reference A can be defined in terms of the vector from frame B, as follows:

The translation and rotation from reference frame {A} to {B}, defined in terms of their respective vectors to point P, can be written in a homogeneous form:

Where is the homogeneous form of vector p of the form

It should be noted:

• This transformation matrix defines the pose of the system, expressing both the translation and rotation with respect to the two frames, such that the negative of the pose is equal to the inverse of the transformation matrix:

• The transformation matrices are non-commutative, therefore, changing the order in which the rotation and translation are performed will result in a different answer altogether.

• Unlike the rotation matrix, the inverse of the transformation matrix is not equal to it’s transpose .

• Homogeneous transformations belong to the special Euclidean group of dimension 2:

Rotation and Translation in 3D

This is not dissimilar to the transformation in 2 dimensions. A vector in 3 dimensions, is defined by a three row vector describing it’s position in x, y and z, e.g.

The vector can also be represented using terms of and in the form of .

The difference between a point and a vector should be kept in mind:

• A point is defined by a vector displacement from another point, typically the origin of a coordinate frame;
• A vector is the difference between two points.

However, both are represented by n real numbers.

Rotation in 3D

For two coordinate frames A and B, above, the transformation in a 3D plane is expressed by:

where the 3×3 matrix is the rotation matrix . In this orthogonal matrix:

• Each column is a unit length vector, orthogonal to all other columns
• As per the 2D case, the inverse of the matrix is the same as the transpose
• The determinant of the matrix is equal to +1, meaning the length of a vector is unchanged by the rotation;
• This matrix belongs to the special orthogonal group of dimension 3:
• Any orientation can be achieved by up to three elemental rotations about x, y or z.

Elementary rotation matrices:

Similarly to the 2D case, Rotation and translation of a coordinate frame in a 3D space can be expressed as a homogeneous transformation. The only difference between the 2D and 3D cases, being the additional elements of rotation and translation about the z-axis which result in the homogeneous transformation being expressed as a 4×4 matrix in the 3 dimensions.

For the transformation between two coordinate frames A and B, with respect to point P, such that

The translation and rotation can be described:

Therefore, the compact form of this can be written:

where is a 3×3 matrix.

Euler Rotation theorem The Swiss mathematician, Leonhard Euler, came up with what’s called the Rotation Theorem, which states that two three-dimensional coordinate frames are related by a sequence of elementary rotations i,e, x, y and z axes. No more than three rotations are needed, but the rotations in this sequence need to be about different axes. Therefore, any rotation can be expressed in terms of no more than three rotations as long as those three rotations are consecutively about different axes.

Euler angles can be used to describe the orientation of a rigid body with respect to a fixed coordinate frame. Proper Euler angles are defined as those angles, where the rotation sequence comprises two of three rotations about the same axes, but not sequentially. For example:

-xyx, zyz, zxz, yxy, xzx and yzy.

Although the sequence zyz is most commonly referred to as Euler angles, all 6 rotations above qualify as Euler angles.

Cardan Angles A sequence of rotations, where each of the three rotations is about a different axes is referred to as a Cardan angle, e.g. xyz, yxz etc.

Gimbal Lock for yaw pitch and roll

Roll, pith and yaw are a common convention used to describe orientation in 3 dimensional space. Roll, pitch and yaw describe the rotation around the x-axis, y-axis and z-axis, respectively.

A Gimbal is a mechanism, commonly used in engineering for navigation systems, which comprises rotatable axis about a fixed point. Gimbals used in navigation systems comprise three axes and are used to maintain a constant frame of reference relative to rotating axes.

In navigation, a gimbal contains what is known as a stabilized platform, which consists of three spinning gyroscopes. that enable the stabilized platform to maintain a fixed orientation with respect to the universe. As the vehicle in which the gimbal is mounted will be be rotating and translating, a series of gimbals are used to connect the stabilized platform to the vehicle in which it is mounted, but in such a way that the motion of the space craft does not introduce any torque or force on the stable platform. This allows the assembly to freely rotate around this axis.

Gimbal lock is a phenomena which occurs when when the pitch reaches radians, where the roll and yaw axis align and it becomes impossible distinguish between the two. This can be represented by a mathematical expression:

, where is parallel to