Optimal design of a Stewart platform using the global transmission index under determinate constraint of workspace
This repository is the application of the Research article first published online October 23, 2017 :
"Optimal design of a Stewart platform using the global transmission index under determinate constraint of workspace " (https://doi.org/10.1177/168781401772088) (https://journals.sagepub.com/doi/10.1177/1687814017720880)
You can use this repository to find the optimal set of parameters which define a Stewart Platform (base_radius,base_angle,platform_radius,platform_angles) based on the maximization of the Global Transmission Index. In this article, screw theory is used to find the transmissability indexes: "The transmissibility indexes describe the transmission ability from the joint force to the platform movement."
The big advantage of this paper is the possibility to optimise the joint configuration of the parallel manipulator based on frame-free and unit homogeneous indexes.
- Due to ne natural frame-dependence of the Jacobian, each index coming from it will be directly dependent on the position in space of the platform : "Thus, those evaluation indexes share common drawbacks such as frame-dependent property and the unit inhomogeneous problem". Because the Transmission Index is based on power transmission, which is a scalar, we can solve such drawbacks: "Thus, this kind of evaluation indexes is frame-free and unit homogeneous."
This repository focuses on replicating the analysis presented in the referenced paper.
First, a brief explanation of screw theory, which serves as the foundation for this paper, is provided. Then, the explained step-by-step procedure outlined in the paper is carefully implemented.
- Explanation of Screw Theory and Global Transmission Index
- Procedure
- Step1: Definition of the Space Docking Mechanism
- Step2: Calculation of Transmission Wrench Screws
- Step3: Obtain the Output Twist Screws of the motion of the platform
- Step4: figure out the Transmission Index when the TWS and OTS have been gotten
- Step5: Fathom out the OATI when the configuration and position of the manipulator are constant
- Step6: work out the GTIc when the configuration of the manipulator is given
- Step7: seek the optimal configuration for manipulator based on the GTIc
Super quick review of Screw Theory.
A screw can be used to describe the movement of a rigid body, usually called twist screw, and also can be used to describe the generalized force applied on a rigid body named wrench screw.
The physical significance of the reciprocal product is the instantaneous power generated by the wrench screw on the rigid body which is moving under the motion described by the twist screw . The reciprocal product between the twist screw St and the wrench screw Sw is the instantaneous power generated by the two screws.
The instantaneous power can be defined as:
For most mechanisms, the prime function is to transmit motion from the input joint to the output joint and the end effector, so the motion of the output joint and the end effector caused by the movement of the input joint represented in screw format is called the output twist screw (OTS). Similarly, the external load exerted on the output joint and the end effector is transmitted to the input joint, so the internal wrench caused by the external load represented in screw format is called the transmission wrench screw (TWS).
The definition of TI (transmission index) for spatial mechanisms is the ratio of instantaneous power and the maximum power generated by the unit TWS (Transmission Wrench Screw Sw) and OTS (Output Twist Screw St) on a rigid body. Hence, the TI between TWS and OTS can be represented as
We can then integrate the transmission index TI , in the whole workspace to obtain a Global Transmission Index relative to the platform joint configuration [r_b, r_p , phi_p , phi_b]. We can then vary the joint configuration to find the configuration which maximises the Global Transmission Index.
In this paper a procedure for finding the optimal configuration for the Space Docking Mechanism is carried out. For the space docking mechanism, the moving platform circle must conform to the international standard, so r_p is set as 0.363m. To calculate the GTI in the whole workspace of the manipulator, we must determine the appropriate zero position first. Here, we also assume the configuration parameter of the manipulator is: [r_b, phi_p , phi_b]= [0.307m, 100°, 60°]. The position of the platform in space is also set to [x, y, z, R, P, Y]= [0, 0.0, 0.2, 0, 0, 0]. Workspace limits are also defined.
# Define parameters
r_p = 0.363; # Radius of platform
r_b = 0.307; # Radius of base
phi_p = 100; # Angle between platform joints
phi_b=60; # Angle between base joints
# Create Stewart Platform instance
platform = StewartPlatform(r_b, phi_b, r_p, phi_p)
pose = [0.0, 0.0, 0.18, 0, 0, 0] # [x, y, z, roll, pitch, yaw] (removed 2 centimeters for platform thickness)
leg_lengths = platform.getIK(pose)
platform.plot()
# Save data points
base_points=platform.b_i
platform_points=platform.p_i
unit_vectors=platform.l_i_unit
workspace_limits = [-0.05,0.05,-0.05,0.05,0.15,0.5] # [x_min, x_max, y_min, y_max, z_min, z_max]
orientation_limits = [-5,5,-5,5,-5,5] # [roll_min, roll_max, pitch_min, pitch_max, yaw_min, yaw_max]
The TWS of the branched chain of the docking mechanism is a pure force whose axis passes through the center of the S pair and U pair (platform leg mechanical joints) simultaneously. And the TWS can be expressed as:
where l_i is the unit vector of the direction of thebranched chain, and r_i is the coordinate of arbitrary point on the branched chain with respect to the reference frame.
# Transmission wrench screw
TWS=np.zeros([6,6])
for i in range(6):
TWS[i,0:3]=unit_vectors[i,:]
TWS[i,3:6]=np.cross(platform_points[i], unit_vectors[i,:])When the Stewart manipulator is only driven by prismatic actuator and other prismatic actuators are locked, the Stewart manipulator becomes a single-DOF manipulator. [...] The five TWSs represented as TWS_j (j=1,...,6 and j=/i), become the constraint wrench of the moving platform, so the OTS_i is reciprocal to all the five TWSs; in this way, one can obtain the following equation:
Here, a kind of method to calculate the reciprocal screw of five given screws is introduced. If you are curious to learn more you can check out the paper itself as the method is a bit complex and long to explain it here. I highly suggest to check it out (:
# Calculating Output Twist Screw (OTS)
OTS=np.zeros([6,6])
TWS_reciprocal=np.zeros([6,6]) # We need the reciprocal of the TWS
A=np.zeros([6,6]) # Matrix for the reciproval of TWS
# Compute the reciprocal of TWS
for i in range(6):
TWS_reciprocal[i,0:3]=TWS[i,3:6]
TWS_reciprocal[i,3:6]=TWS[i,0:3]
A=np.copy(TWS_reciprocal) # Matrix for the reciproval of TWS
# Algorithm for calculating OTS
for k in range(5, -1, -1):
A_without_row=np.copy(np.delete(A, k, axis = 0)) # Each k iteration remove a row from the matrix.
det_vect=np.zeros(6) # Determinant vector
for i in range(6):
A_local=np.copy(A_without_row)
A_local = np.delete(A_local, i, axis = 1) # Each i iteration remove i column
det_vect[i]=np.linalg.det(A_local) # Compute determinant of remaining 5x5 matrix
gain=1/np.linalg.norm(np.array([-det_vect[0],det_vect[1],-det_vect[2]])) # Calculate gain
OTS[k,:]=gain*np.array([-det_vect[0],det_vect[1],-det_vect[2],det_vect[3],-det_vect[4],det_vect[5]]) # Find OTS vector for each k iterationYou can check if the method has worked correctly by verifying the reciprocity of the two vectors. If the two vectors are reciprocal it means that their dot product will give 0 as a result.
#TWSi=[w,v] this lays in each joint (force from legs)
#OTSi=[w,v] this lays in the platform point center (velocity from legs on the P point)
print("Check if OTS and TWS are reciprocal with j=/i \n")
# if j =/ i
j=0;
i=1;
print("Dot product of reciprocal vectors is 0")
print("np.dot(TWS_reciprocal[j],OTS[i]): ",np.dot(TWS_reciprocal[j],OTS[i]),"\n")
print("Dot product of non-reciprocal vectors is different from 0")
print("np.dot(TWS_reciprocal[j],OTS[j]): ",np.dot(TWS_reciprocal[j],OTS[j]))Check if OTS and TWS are reciprocal with j=/i
Dot product of reciprocal vectors is 0
np.dot(TWS_reciprocal[j],OTS[i]): 5.551115123125783e-17
Dot product of non-reciprocal vectors is different from 0
np.dot(TWS_reciprocal[j],OTS[j]): 0.5920283360375702
To calculate the Transmission Index TI we need to find the numerator and the denominator of such equation:
The numerator can be easily found with the dot product of TWS_j with OTS_j, for the denominator we can look directly for the max:
When the two screws , TWS and OTS , are given, their pitch, h1 and h2,are thought to be constant.
To calculate d_max we firstly have to find the distance from the characteristic point of the TWS (joint platform attachment) to the axis line of the OTS. To do so we need to find the equation of the points laying on OTS vector.
#### Calculating Transmission Index and returning the minimum
w_OTS=np.copy(OTS[:,0:3]) # Unit vector of the rotational velocity of OTS
v_OTS=np.copy(OTS[:,3:6]) # The dual of the unit screw of OTS
d_max=np.zeros(6) # Stores max distance
d=np.zeros([6,3]) # " where d is the vector from the characteristic point Pi to arbitrary point on the axis line of the OTSi"
h1_vect=np.zeros(6) # pitch of the unit screw of TWS
h2_vect=np.zeros(6) # pitch of the unit screw of OTS
TI_vect=np.ones([6])*100 # Initializing Trasmission Index Vector
# We need a point Laying on each OTS vector.
# r x w_i = v_i - h*w_i (screw theory equation defining OTS)
# h = v_i * w_i
# r x w_i = v_i - (v_i * w_i) * w_i
r_OTS = np.zeros([6,3]) # Points will contain the point laying on OTS vector
c = np.zeros([6,3]) # Right side of cross product
h = np.zeros(6) # Pitch of OTS[i]
t=0; # Line parameter, value can be chosen freely (we will do cross products, the point has to just stay on the OTS vector)
# Calculate TI for each joint
for i in range(6):
# Calculate points laying on OTSi vector
# r x w_i = v_i - h*w_i
# h = v_i * w_i
# r x w_i = v_i - (v_i * w_i) * w_i
# a = np.cross(b,c)/np.dot(b,b)+t*b is a solution for all t. if a and c are orthogonal
# r_i = np.cross(w_OTS[i],c[i])/np.dot(w_OTS[i],w_OTS[i])+t*w_OTS[i] , t=0 --> r_i = np.cross(w_OTS[i],c[i])
# where c[i] = v_i - (v_i * w_i) * w_i
h[i]=np.dot(w_OTS[i],v_OTS[i]) # pitch of OTS[i]
c[i]= v_OTS[i] - h[i]*w_OTS[i] # right side of cross product r x w_i = v_i - (v_i * w_i) * w_i
r_OTS[i]=np.cross(w_OTS[i],c[i]) +t*w_OTS[i] # points laying on OTSi
# Calculate max perpendicular distance from each platform point to their respective OTS vector
d[i] = platform_points[i]-r_OTS[i] # " [...] where d is the vector from the characteristic point Pi to arbitrary point on the axis line of the OTSi"
d_max[i]=np.linalg.norm(np.cross(d[i],w_OTS[i])) # max perpendicular distance from Pi to the OTS vector
# Calculate unit pitch of TWSi and OTSi
h1_vect[i]=np.dot(TWS_reciprocal[i,0:3],TWS_reciprocal[i,3:6]) # unit pitch of TWSi
h2_vect[i]=np.dot(OTS[i,3:6],OTS[i,0:3]) # unit pitch of OTSi
# Calculating TI for each joint
numerator=np.abs(np.dot(TWS_reciprocal[i],OTS[i]))
denominator=np.sqrt((h1_vect[i] + h2_vect[i]) ** 2 + d_max[i] ** 2)
TI_vect[i]=numerator/denominator
# Taking the minimum Transmission Index value.
TI_min=np.min(TI_vect)As we need to plot TI for every orientation in our workspace, we will create a function getTI which incorporates step2, step3 and step4. We can then calculate the orientation distribution of TI (OATI Orientation Avarage Transmission Index) in the whole discretized orientation workspace.
position=pose[0:3] # pose from the definition of the Space Docking Mechanism
# Define discretization space
N=20 # Discretization
roll_min, roll_max, pitch_min, pitch_max, yaw_min, yaw_max = orientation_limits # Extracting orientation_limits
roll_vect = np.linspace(roll_min, roll_max, N)
pitch_vect = np.linspace(pitch_min, pitch_max, N)
yaw_vect = np.linspace(yaw_min, yaw_max, N)
rr, pp, yy = np.meshgrid(roll_vect, pitch_vect, yaw_vect, indexing='ij')
orientations = np.vstack([rr.ravel(), pp.ravel(), yy.ravel()]).T
# orientation average transmission index (OATI)
Holder = []
for orient in orientations:
p_vect = np.hstack((position, orient))
platform.getIK(p_vect)
TI=getTI(platform)
Holder.append(np.append(orient, TI))
Holder=np.array(Holder)
We can then calculate the OATI distribution in the position workspace
N= 4 # Define discretization of space
# Define discretization space for positions
x_min, x_max, y_min, y_max, z_min, z_max = workspace_limits
N_pos = N
x_vect = np.linspace(x_min, x_max, N_pos)
y_vect = np.linspace(y_min, y_max, N_pos)
z_vect = np.linspace(z_min, z_max, N_pos)
xx, yy, zz = np.meshgrid(x_vect, y_vect, z_vect, indexing='ij')
positions = np.vstack([xx.ravel(), yy.ravel(), zz.ravel()]).T
# Define discretization space for orientations
roll_min, roll_max, pitch_min, pitch_max, yaw_min, yaw_max = orientation_limits
N_orient = N
roll_vect = np.linspace(roll_min, roll_max, N_orient)
pitch_vect = np.linspace(pitch_min, pitch_max, N_orient)
yaw_vect = np.linspace(yaw_min, yaw_max, N_orient)
rr, pp, yy = np.meshgrid(roll_vect, pitch_vect, yaw_vect, indexing='ij')
orientations = np.vstack([rr.ravel(), pp.ravel(), yy.ravel()]).T
# Iterate through each position and each orientation
Holder_mean_orient = []
Holder = []
for pos in positions:
Holder_TI = []
for orient in orientations:
p_vect = np.hstack((pos, orient))
platform.getIK(p_vect)
TI=getTI(platform)
Holder_TI.append(getTI(platform))# get every TI in orientations
OATI=np.mean(Holder_TI)
Holder_mean_orient.append(OATI) #AOTI
Holder.append(np.append(pos, OATI))
GTI=np.mean(Holder_mean_orient)
Holder=np.array(Holder)
We can write a function that takes calculates the Global Transmission Index for every try of platform variables.
# define platform parameters
phi_p_min, phi_p_max, phi_b_min, phi_b_max, r_b_min, r_b_max = [10,110,10,110,0.3,0.45] # define platform variables min and max
N_plat = 6 # N x N x N # define number of platform to test N_plat^3
phi_p_vect = np.linspace(phi_p_min, phi_p_max, N_plat)
phi_b_vect = np.linspace(phi_b_min, phi_b_max, N_plat)
r_b_vect = np.linspace(r_b_min, r_b_max, N_plat)
N_WS = 4 # discretize workspaces
xx, yy, zz = np.meshgrid(phi_p_vect, phi_b_vect, r_b_vect, indexing='ij')
platforms = np.vstack([xx.ravel(), yy.ravel(), zz.ravel()]).T
tot=len(platforms)
Holder = []
i=1
# Iterate through each platform
for plat in platforms:
phi_p = plat[0]; # select platform angle
phi_b = plat[1];# select base angle
r_b = plat[2]; # select radius of base
r_p = 0.363; # fixed by constraints
platform = StewartPlatform(r_b, phi_b, r_p, phi_p) # create platform object
GTI = getGTI(platform,workspace_limits,orientation_limits,N_WS) # calculate GTI for platform object
Holder.append(np.append(plat, GTI))
print("number of platform tested:", i ,"/", tot)
i+=1
Holder=np.array(Holder)
We can notice the instability of the configuration when the base angle and the platform angle are equal (the platform is inherently unstable).
We can also plot the GTI surfaces to better see the change.
We can then select the platform with the highest Global Transmission Index
optimal = np.max(Holder[:,3]) # find max GTI
config=np.where(Holder[:,3]==optimal)
optimal_config=Holder[config[0][0],:] # select platform with max GTI
phi_p = optimal_config[0]
phi_b = optimal_config[1]
r_b = optimal_config[2]
r_p = 0.363
N_WS = 5
optimal_platform=StewartPlatform(r_b,phi_b,r_p,phi_p)
optimal_platform.getIK([0,0,0.3,0,0,0])
optimal_platform.plot()
print("optimal phi_p", optimal_config[0])
print("optimal phi_b", optimal_config[1])
print("optimal r_b", optimal_config[2])
print("GTI: ",optimal_config[3])
optimal phi_p 30.0
optimal phi_b 110.0
optimal r_b 0.45
GTI: 0.900375823640452