使用Phantom omni力反馈设备控制机器人

传统的工业机器人普遍采用电机 、齿轮减速器 、关节轴三者直接连接的传动机构，这种机构要求电机与减速器安装在机械臂关节附近，其缺点是对于多关节机械臂，下一级关节的电机与减速器等驱动装置成为上一级关节的额外负载 。这一额外负载带来的负面影响往往超过机械臂连杆等必要结构件，因此提高了对机械臂动力和驱动元件的要求，由此造成整体重量 、体积 、造价和内部消耗的增加，降低了机械臂对外做功的能力和效率。为了避免这一问题许多主动式医疗机器人采用绳驱动方式，使用柔性绳索远距离的传递运动和力矩。将驱动装置安装到远离相关关节的基座上，使得整个机械臂的结构紧凑 ，在减轻机械臂整体重量的同时提高了对外做功的能力。下图所示的就是达芬奇手术机器人的末端夹子，采用绳驱动方式可以在较小的机械结构内实现多自由度的灵活运动：

下面尝试在VREP中使用力反馈设备Phantom omni来控制医疗机器人的末端夹子。机构简图如下图所示，一共有5个自由度：其中移动自由度负责进给，一个整体旋转自由度，还有两个左右、上下弯曲的自由度，最后是控制夹子张合的自由度。

删除Solidworks CAD模型中的一些不必要特征，比如倒角、内部孔、螺纹孔等，导出成STL文件，再导入到VREP中。然后还需要进行继续化简，减少网格数量，当网格数减少到不影响几何外观就可以了。接下来在相应的位置添加关节，设置关节运动范围（软件限位），并将其设为Inverse kinematics模式。

设置Calculation Modules中的Inverse kinematics：

搭建好模型后可以先用键盘进行测试，按方向键移动ikTarget，VREP会根据构型自动计算运动学逆解，然后将ikTip移动到ikTarget处（这样就会带着整个机构运动）。

键盘控制的脚本代码如下：

• if (sim_call_type==sim_childscriptcall_initialization) then
• -- Put some initialization code here
• targetHandle = simGetObjectHandle('ikTarget')
• end
• if (sim_call_type==sim_childscriptcall_actuation) then
• -- Put your main ACTUATION code here
• end
• if (sim_call_type==sim_childscriptcall_sensing) then
• -- Put your main SENSING code here
• -- Read the keyboard messages (make sure the focus is on the main window, scene view):
• message, auxiliaryData = simGetSimulatorMessage()
• if (message == sim_message_keypress) then
• if (auxiliaryData[1]==119) then
• -- W key
• local p = simGetObjectPosition(targetHandle, -1)
• p[1] = p[1] - 0.001
• simSetObjectPosition(targetHandle, -1, p)
• end
• if (auxiliaryData[1]==115) then
• -- S key
• local p = simGetObjectPosition(targetHandle, -1)
• p[1] = p[1] + 0.001
• simSetObjectPosition(targetHandle, -1, p)
• end
• if (auxiliaryData[1]==2007) then
• -- up key
• local p = simGetObjectPosition(targetHandle, -1)
• p[3] = p[3] + 0.001
• simSetObjectPosition(targetHandle, -1, p)
• end
• if (auxiliaryData[1]==2008) then
• -- down key
• local p = simGetObjectPosition(targetHandle, -1)
• p[3] = p[3] - 0.001
• simSetObjectPosition(targetHandle, -1, p)
• end
• if (auxiliaryData[1]==2009) then
• -- left key
• local p = simGetObjectPosition(targetHandle, -1)
• p[2] = p[2] - 0.001
• simSetObjectPosition(targetHandle, -1, p)
• end
• if (auxiliaryData[1]==2010) then
• -- right key
• local p = simGetObjectPosition(targetHandle, -1)
• p[2] = p[2] + 0.001
• simSetObjectPosition(targetHandle, -1, p)
• end
• end
• end
• if (sim_call_type==sim_childscriptcall_cleanup) then
• -- Put some restoration code here
• end

测试没问题后可以使用CHAI3D插件来连接力反馈设备。这里采用增量控制的模式，即计算当前时刻与前一时刻手柄位置在X、Y、Z方向上的差，然后控制VREP中的ikTarget控制点按相应的增量移动。注意在VREP中机器人向前运动是沿X轴负方向、向上运动是沿Z轴正方向、向右运动是沿Y轴正方向，这与CHAI3D中坐标系的定义一致（向前推手柄是沿着X轴负方向...），因此可以使用力反馈设备直观的控制机器人的运动。当然如果坐标系定义不一致，需要进行简单的转换才行。

下面的代码在sensing探测部分会使用simExtCHAI3D_readPosition来读取当前操作手柄的位置和按钮状态。按照VREP默认设置，这部分代码会50ms执行一次，这里会出现一个问题：如果采样速率太快，会导致前后两次采集到的位置数据偏差为零（人手的操作频率没那么快，还来不及改变位置），那么输出的控制量就一直是零，这样就没办法用增量控制的方式来操控机器人。解决办法是延迟几个周期再采样，等到有足够大的偏差之后再生成控制量。还有一个问题是使用CHAI3D返回的数据以“米”为单位，而VREP世界中的单位有可能未知，那么使用增量控制时需要对控制量乘一个比例系数，避免因操作端微小的移动造成从动端运动量过大，超出关节限制（无法到达的逆解）。或者可以调节比例系数，用操作端的大位移来控制从动端的小位移，实现精细控制。

• if (sim_call_type==sim_childscriptcall_initialization) then
• -- Check if the plugin is loaded:
• moduleName=0
• moduleVersion=0
• index=0
• pluginNotFound=true
• while moduleName do
• moduleName,moduleVersion=simGetModuleName(index)
• if (moduleName=='CHAI3D') then
• pluginNotFound=false
• end
• index=index+1
• end
• if (pluginNotFound) then
• simDisplayDialog('Error','CHAI3D plugin was not found, or was not correctly initialized (v_repExtCHAI3D.dll).',sim_dlgstyle_ok,false,nil,{0.8,0,0,0,0,0},{0.5,0,0,1,1,1})
• else
• -- Start the device:
• local toolRadius = 0.001 -- the radius of the tool
• local workspaceRadius = 0.2 -- the workspace radius
• if simExtCHAI3D_start(0, toolRadius,workspaceRadius) ~= 1 then
• simDisplayDialog('Error','Device failed to initialize.',sim_dlgstyle_ok,false,nil,{0.8,0,0,0,0,0},{0.5,0,0,1,1,1})
• else
• CHAI3DPluginInitialized = true
• end
• end
• targetHandle = simGetObjectHandle('ikTarget')
• deltaPos = {0, 0, 0}
• counter = 0
• ratio = 50
• end
• if (sim_call_type==sim_childscriptcall_actuation) then
• if buttonState ==1 then -- press the button
• local p = simGetObjectPosition(targetHandle, -1) -- get the target position
• -- add increment of the tool tip
• p[1] = p[1] + deltaPos[1] / ratio
• p[2] = p[2] + deltaPos[2] / ratio
• p[3] = p[3] + deltaPos[3] / ratio
• simSetObjectPosition(targetHandle, -1, p) -- move to the absolute position
• end
• end
• if (sim_call_type==sim_childscriptcall_sensing) then
• if CHAI3DPluginInitialized then
• -- Read the current position of the cursor:
• local currentToolPosition = simExtCHAI3D_readPosition(0)
• -- Read the buttons of the device:
• buttonState = simExtCHAI3D_readButtons(0)
• counter = counter + 1 -- increase the counter
• if counter % 30 == 1 then -- keep the value
• prevToolPosition = currentToolPosition
• end
• if counter % 30 == 0 then -- calculate tool tip increment
• deltaPos[1] = currentToolPosition[1] - prevToolPosition[1] -- X-axis increment
• deltaPos[2] = currentToolPosition[2] - prevToolPosition[2] -- Y-axis increment
• deltaPos[3] = currentToolPosition[3] - prevToolPosition[3] -- Z-axis increment
• counter = 0 -- reset counter
• local info = string.format("CurrentPosition:%.2f,%.2f,%.2f DeltaPosition:%.2f,%.2f,%.2f",
• currentToolPosition[1],currentToolPosition[2],currentToolPosition[3],
• deltaPos[1],deltaPos[2],deltaPos[3])
• simAddStatusbarMessage(info)
• end
• end
• end
• if (sim_call_type==sim_childscriptcall_cleanup) then
• if CHAI3DPluginInitialized then
• -- Disconnects all devices and removes all objects from the scene
• simExtCHAI3D_reset()
• end
• end

将力反馈设备手柄移动到合适的位置之后就可以按住按钮开始操控机器人，松开按钮会停止控制。如果在VREP虚拟场景中添加其它物体（比如障碍物），则还可以模拟环境中的力（接触力、重力、摩擦力、弹簧力等）让操控着“感觉”到。如果实际机器人上装有力传感器，则在用Phantom omni控制机器人的同时也能读取力的信息，反馈给操作者。

下面是使用Phantom omni来控制机器人的动态图，黄色的轨迹为使用Graph记录的控制点的空间位置：

对于该机构也可以自己实现运动学逆解的数值算法，下面给出伪逆矩阵法和阻尼最小二乘法的参考：

• import math
• import numpy as np
• # link length
• L1 = 1
• L2 = 1
• gamma = 1 # step size
• lamda = 0.2 # damping constant (DLS-method)
• stol = 1e-3 # tolerance
• nm = 10 # initial error
• count = 0 # iteration count
• ilimit = 20 # maximum iteration
• # numerical method for inverse kinematics
• method = 'Pseudo Inverse' # 'Pseudo Inverse', 'DLS', ...
• # initial joint value
• q = np.array([0, 0, math.pi/2, 0]) # [theta1, d1, theta2, theta3]
• # target position
• target_pos = np.array([1, 0, 2]) # [x,y,z]
• while True:
• if(nm > stol):
• # forward kinematics:
• x = np.array([math.cos(q[0])*math.cos(q[2])*(L1+L2*math.cos(q[3]))+L2*math.sin(q[0])*math.sin(q[3]),
• math.cos(q[2])*math.sin(q[0])*(L1+L2*math.cos(q[3]))-L2*math.cos(q[0])*math.sin(q[3]),
• q[1]+(L1+L2*math.cos(q[3]))*math.sin(q[2])])
• # compute error
• error = target_pos - x
• # compute Jacobian
• J11 = -math.sin(q[0])*math.cos(q[2])*(L1+L2*math.cos(q[3]))+L2*math.cos(q[0])*math.sin(q[3])
• J12 = 0
• J13 = -math.sin(q[2])*math.cos(q[0])*(L1+L2*math.cos(q[3]))
• J14 = L2*(math.sin(q[0])*math.cos(q[3])-math.cos(q[0])*math.cos(q[2])*math.sin(q[3]))
• J21 = math.cos(q[0])*math.cos(q[2])*(L1+L2*math.cos(q[3]))+L2*math.sin(q[0])*math.sin(q[3])
• J22 = 0
• J23 = -math.sin(q[0])*math.sin(q[2])*(L1+L2*math.cos(q[3]))
• J24 = -L2*(math.cos(q[0])*math.cos(q[3])+math.sin(q[0])*math.cos(q[2])*math.sin(q[3]))
• J31 = 0
• J32 = 1
• J33 = math.cos(q[2])*(L1+L2*math.cos(q[3]))
• J34 = -L2*math.sin(q[2])*math.sin(q[3])
• J = np.array([[J11,J12,J13,J14],[J21,J22,J23,J24],[J31,J32,J33,J34]])
• if method == 'Pseudo Inverse':
• # Pseudo Inverse Method
• J_pseudo = np.dot(J.transpose(), np.linalg.inv(J.dot(J.transpose()))) # compute pseudo inverse
• dq = gamma * J_pseudo.dot(error)
• if method == 'DLS':
• # Damped Least Squares Method
• f = np.linalg.solve(J.dot(J.transpose())+lamda**2*np.identity(3), error)
• dq = np.dot(J.transpose(), f)
• # update joint position
• q = q + dq
• nm = np.linalg.norm(error)
• count = count + 1
• if count > ilimit:
• print "Solution wouldn't converge!"
• break
• else:
• # print result
• print 'theta1=' + str(q[0]*180/math.pi)+' d1='+str(q[1]) + ' theta2=' + str((q[2]-math.pi/2)*180/math.pi)+' theta3=' + str(q[3]*180/math.pi)
• print 'Current position: %.2f, %.2f, %.2f' %(x[0],x[1],x[2])
• print str(count)+' iterations'+' err:'+str(nm)
• break