Code:

import pybullet as pimport timeimport numpy as npfrom numpy import linalg as laimport mathp.connect(p.GUI)p.createCollisionShape(p.GEOM_PLANE)p.createMultiBody(0,0)# PARAMETERSn_links = 30          # num linksdx_link = 0.02        # length of link segmentlink_mass = 0.005     # 5gbase_mass = 0.1       # 100gjoint_friction = 0.0005  # rotational joint friction [N/(rad/s)]# setup shapeslink_shape = p.createCollisionShape(p.GEOM_CYLINDER,                                     radius=0.005,                                     height=dx_link,                                     collisionFramePosition=[0, 0, -dx_link/2])base_shape = p.createCollisionShape(p.GEOM_BOX,                                    halfExtents=[0.01, 0.01, 0.01])linkMasses = [link_mass]*n_linkslinkCollisionShapeIndices=[link_shape]*n_linkslinkVisualShapeIndices=[-1]*n_links# relative position linkPositions = []for i in range(n_links):    linkPositions.append([0,0,-dx_link])# cm positionslinkOrientations=[[0,0,0,1]]*n_linkslinkInertialFramePositions=[[0,0,0]]*n_linkslinkInertialFrameOrientations=[[0,0,0,1]]*n_links# connection graphindices = range(n_links)# use spherical jointsjointTypes = [p.JOINT_SPHERICAL]*n_linksjointTypes[1] = p.JOINT_FIXED# rotational axis (dosnt't matter, spherical)axis = [[1,0,0]]*n_links# create rope bodyvisualShapeId = -1basePosition = [0,0,2]baseOrientation = [0,0,0,1]rope = p.createMultiBody(base_mass,base_shape,visualShapeId,basePosition,baseOrientation,                        linkMasses=linkMasses,                        linkCollisionShapeIndices=linkCollisionShapeIndices,                        linkVisualShapeIndices=linkVisualShapeIndices,                        linkPositions=linkPositions,                        linkOrientations=linkOrientations,                        linkInertialFramePositions=linkInertialFramePositions,                        linkInertialFrameOrientations=linkInertialFrameOrientations,                        linkParentIndices=indices,                        linkJointTypes=jointTypes,                        linkJointAxis=axis)                        #flags=p.URDF_USE_SELF_COLLISION)n_joints = p.getNumJoints(rope)# remove stiffness in motors, add friction forcefriction_vec = [joint_friction]*3   # same all axiscontrol_mode = p.POSITION_CONTROL   # set pos control modefor j in range(n_joints):  p.setJointMotorControlMultiDof(rope,j,control_mode,                                 targetPosition=[0,0,0,1],                                 targetVelocity=[0,0,0],                                 positionGain=0,                                 velocityGain=1,                                 force=friction_vec)p.setGravity(0,0,-9.81)# fixed constrain to keep root cube in placeroot_robe_c = p.createConstraint(rope, -1, -1, -1, p.JOINT_FIXED, [0, 0, 0], [0, 0, 0], [0, 0, 2])# some traj to inject motionamplitude_x = 0.3amplitude_y = 0.0freq = 0.6# manually simulate joint dampingDamping = 0.001t = 0.0freq_sim = 240while 1:  time.sleep(1./freq_sim)  t += 1./freq_sim  # some trajectory  ux = amplitude_x*math.sin(2*math.pi*freq*t)  uy = amplitude_y*math.cos(2*math.pi*freq*t)    # move base arround  pivot = [ux, uy, 2]  orn = p.getQuaternionFromEuler([0, 0, 0])  p.changeConstraint(root_robe_c, pivot, jointChildFrameOrientation=orn, maxForce=500)  # manually apply viscous friction: f_damping = -Damping*omega  '''  for j in range(n_joints):    q_state = p.getJointStateMultiDof(rope, j)    omega = q_state[1]    f_dampling = [-Damping*v for v in omega]    p.setJointMotorControlMultiDof(rope, j, p.TORQUE_CONTROL, force=f_dampling)  '''  # update  p.stepSimulation()