|
8 | 8 | import time |
9 | 9 | import subprocess |
10 | 10 | import sys |
11 | | -sys.path += [os.getcwd()+'/../exercizes'] |
| 11 | + |
| 12 | +sys.path += [os.getcwd() + "/../exercizes"] |
12 | 13 | import plot_utils as plut |
13 | 14 | import matplotlib.pyplot as plt |
14 | 15 |
|
15 | 16 | np.set_printoptions(precision=3, linewidth=200, suppress=True) |
16 | 17 |
|
17 | 18 | LINE_WIDTH = 60 |
18 | | -print("".center(LINE_WIDTH,'#')) |
19 | | -print(" Test TSID with Quadruped Robot ".center(LINE_WIDTH, '#')) |
20 | | -print("".center(LINE_WIDTH,'#'), '\n') |
21 | | - |
22 | | -mu = 0.3 # friction coefficient |
23 | | -fMin = 1.0 # minimum normal force |
24 | | -fMax = 100.0 # maximum normal force |
25 | | -contact_frames = ['BL_contact', 'BR_contact', 'FL_contact', 'FR_contact'] |
26 | | -contactNormal = np.array([0., 0., 1.]) # direction of the normal to the contact surface |
27 | | - |
28 | | -w_com = 1.0 # weight of center of mass task |
29 | | -w_posture = 1e-3 # weight of joint posture task |
30 | | -w_forceRef = 1e-5 # weight of force regularization task |
31 | | - |
32 | | -kp_contact = 10.0 # proportional gain of contact constraint |
33 | | -kp_com = 10.0 # proportional gain of center of mass task |
34 | | -kp_posture = 10.0 # proportional gain of joint posture task |
35 | | - |
36 | | -dt = 0.001 # controller time step |
37 | | -PRINT_N = 500 # print every PRINT_N time steps |
38 | | -DISPLAY_N = 25 # update robot configuration in viwewer every DISPLAY_N time steps |
39 | | -N_SIMULATION = 6000 # number of time steps simulated |
| 19 | +print("".center(LINE_WIDTH, "#")) |
| 20 | +print(" Test TSID with Quadruped Robot ".center(LINE_WIDTH, "#")) |
| 21 | +print("".center(LINE_WIDTH, "#"), "\n") |
| 22 | + |
| 23 | +mu = 0.3 # friction coefficient |
| 24 | +fMin = 1.0 # minimum normal force |
| 25 | +fMax = 100.0 # maximum normal force |
| 26 | +contact_frames = ["BL_contact", "BR_contact", "FL_contact", "FR_contact"] |
| 27 | +contactNormal = np.array( |
| 28 | + [0.0, 0.0, 1.0] |
| 29 | +) # direction of the normal to the contact surface |
| 30 | + |
| 31 | +w_com = 1.0 # weight of center of mass task |
| 32 | +w_posture = 1e-3 # weight of joint posture task |
| 33 | +w_forceRef = 1e-5 # weight of force regularization task |
| 34 | + |
| 35 | +kp_contact = 10.0 # proportional gain of contact constraint |
| 36 | +kp_com = 10.0 # proportional gain of center of mass task |
| 37 | +kp_posture = 10.0 # proportional gain of joint posture task |
| 38 | + |
| 39 | +dt = 0.001 # controller time step |
| 40 | +PRINT_N = 500 # print every PRINT_N time steps |
| 41 | +DISPLAY_N = 25 # update robot configuration in viwewer every DISPLAY_N time steps |
| 42 | +N_SIMULATION = 6000 # number of time steps simulated |
40 | 43 |
|
41 | 44 | filename = str(os.path.dirname(os.path.abspath(__file__))) |
42 | | -path = filename + '/../models' |
43 | | -urdf = path + '/quadruped/urdf/quadruped.urdf' |
| 45 | +path = filename + "/../models" |
| 46 | +urdf = path + "/quadruped/urdf/quadruped.urdf" |
44 | 47 | vector = pin.StdVec_StdString() |
45 | 48 | vector.extend(item for item in path) |
46 | 49 | robot = tsid.RobotWrapper(urdf, vector, pin.JointModelFreeFlyer(), False) |
47 | 50 |
|
48 | 51 | # for gepetto viewer |
49 | | -robot_display = pin.RobotWrapper.BuildFromURDF(urdf, [path, ], pin.JointModelFreeFlyer()) |
| 52 | +robot_display = pin.RobotWrapper.BuildFromURDF( |
| 53 | + urdf, |
| 54 | + [ |
| 55 | + path, |
| 56 | + ], |
| 57 | + pin.JointModelFreeFlyer(), |
| 58 | +) |
50 | 59 | l = subprocess.getstatusoutput("ps aux |grep 'gepetto-gui'|grep -v 'grep'|wc -l") |
51 | 60 | if int(l[1]) == 0: |
52 | | - os.system('gepetto-gui &') |
| 61 | + os.system("gepetto-gui &") |
53 | 62 | time.sleep(1) |
54 | 63 | cl = gepetto.corbaserver.Client() |
55 | 64 | gui = cl.gui |
56 | 65 | robot_display.initViewer(loadModel=True) |
57 | 66 |
|
58 | | -#model = robot.model() |
59 | | -#pin.loadReferenceConfigurations(model, srdf, False) |
60 | | -#q = model.referenceConfigurations["half_sitting"] |
61 | | -#q = pin.getNeutralConfigurationFromSrdf(robot.model(), srdf, False) |
62 | | -#q = robot_display.model.neutralConfiguration #np.zeros(robot.nq) |
| 67 | +# model = robot.model() |
| 68 | +# pin.loadReferenceConfigurations(model, srdf, False) |
| 69 | +# q = model.referenceConfigurations["half_sitting"] |
| 70 | +# q = pin.getNeutralConfigurationFromSrdf(robot.model(), srdf, False) |
| 71 | +# q = robot_display.model.neutralConfiguration #np.zeros(robot.nq) |
63 | 72 | q = np.zeros(robot.nq) |
64 | 73 | q[6] = 1.0 |
65 | 74 | q[2] += 0.5 |
66 | 75 | for i in range(4): |
67 | | - q[7 + 2*i] = -0.8 |
68 | | - q[8 + 2*i] = 1.6 |
| 76 | + q[7 + 2 * i] = -0.8 |
| 77 | + q[8 + 2 * i] = 1.6 |
69 | 78 | v = np.zeros(robot.nv) |
70 | 79 |
|
71 | 80 | robot_display.displayCollisions(False) |
|
74 | 83 |
|
75 | 84 | assert [robot.model().existFrame(name) for name in contact_frames] |
76 | 85 |
|
77 | | -t = 0.0 # time |
| 86 | +t = 0.0 # time |
78 | 87 | invdyn = tsid.InverseDynamicsFormulationAccForce("tsid", robot, False) |
79 | 88 | invdyn.computeProblemData(t, q, v) |
80 | 89 | data = invdyn.data() |
|
84 | 93 | q[2] -= robot.framePosition(data, id_contact).translation[2] |
85 | 94 | robot.computeAllTerms(data, q, v) |
86 | 95 |
|
87 | | -contacts = 4*[None] |
| 96 | +contacts = 4 * [None] |
88 | 97 | for i, name in enumerate(contact_frames): |
89 | | - contacts[i] =tsid.ContactPoint(name, robot, name, contactNormal, mu, fMin, fMax) |
| 98 | + contacts[i] = tsid.ContactPoint(name, robot, name, contactNormal, mu, fMin, fMax) |
90 | 99 | contacts[i].setKp(kp_contact * np.ones(3)) |
91 | 100 | contacts[i].setKd(2.0 * np.sqrt(kp_contact) * np.ones(3)) |
92 | 101 | H_rf_ref = robot.framePosition(data, robot.model().getFrameId(name)) |
|
100 | 109 | invdyn.addMotionTask(comTask, w_com, 1, 0.0) |
101 | 110 |
|
102 | 111 | postureTask = tsid.TaskJointPosture("task-posture", robot) |
103 | | -postureTask.setKp(kp_posture * np.ones(robot.nv-6)) |
104 | | -postureTask.setKd(2.0 * np.sqrt(kp_posture) * np.ones(robot.nv-6)) |
| 112 | +postureTask.setKp(kp_posture * np.ones(robot.nv - 6)) |
| 113 | +postureTask.setKd(2.0 * np.sqrt(kp_posture) * np.ones(robot.nv - 6)) |
105 | 114 | invdyn.addMotionTask(postureTask, w_posture, 1, 0.0) |
106 | 115 |
|
107 | 116 | com_ref = robot.com(data) |
|
111 | 120 | q_ref = q[7:] |
112 | 121 | trajPosture = tsid.TrajectoryEuclidianConstant("traj_joint", q_ref) |
113 | 122 |
|
114 | | -print("Create QP solver with ", invdyn.nVar, " variables, ", invdyn.nEq, |
115 | | - " equality and ", invdyn.nIn, " inequality constraints") |
| 123 | +print( |
| 124 | + "Create QP solver with ", |
| 125 | + invdyn.nVar, |
| 126 | + " variables, ", |
| 127 | + invdyn.nEq, |
| 128 | + " equality and ", |
| 129 | + invdyn.nIn, |
| 130 | + " inequality constraints", |
| 131 | +) |
116 | 132 | solver = tsid.SolverHQuadProgFast("qp solver") |
117 | 133 | solver.resize(invdyn.nVar, invdyn.nEq, invdyn.nIn) |
118 | 134 |
|
119 | | -com_pos = np.empty((3, N_SIMULATION))*nan |
120 | | -com_vel = np.empty((3, N_SIMULATION))*nan |
121 | | -com_acc = np.empty((3, N_SIMULATION))*nan |
122 | | - |
123 | | -com_pos_ref = np.empty((3, N_SIMULATION))*nan |
124 | | -com_vel_ref = np.empty((3, N_SIMULATION))*nan |
125 | | -com_acc_ref = np.empty((3, N_SIMULATION))*nan |
126 | | -com_acc_des = np.empty((3, N_SIMULATION))*nan # acc_des = acc_ref - Kp*pos_err - Kd*vel_err |
127 | | - |
128 | | -offset = robot.com(data) + np.array([0.0, 0.0, 0.0]) |
129 | | -amp = np.array([0.0, 0.03, 0.0]) |
130 | | -two_pi_f = 2*np.pi*np.array([0.0, 2.0, 0.7]) |
131 | | -two_pi_f_amp = np.multiply(two_pi_f,amp) |
| 135 | +com_pos = np.empty((3, N_SIMULATION)) * nan |
| 136 | +com_vel = np.empty((3, N_SIMULATION)) * nan |
| 137 | +com_acc = np.empty((3, N_SIMULATION)) * nan |
| 138 | + |
| 139 | +com_pos_ref = np.empty((3, N_SIMULATION)) * nan |
| 140 | +com_vel_ref = np.empty((3, N_SIMULATION)) * nan |
| 141 | +com_acc_ref = np.empty((3, N_SIMULATION)) * nan |
| 142 | +com_acc_des = ( |
| 143 | + np.empty((3, N_SIMULATION)) * nan |
| 144 | +) # acc_des = acc_ref - Kp*pos_err - Kd*vel_err |
| 145 | + |
| 146 | +offset = robot.com(data) + np.array([0.0, 0.0, 0.0]) |
| 147 | +amp = np.array([0.0, 0.03, 0.0]) |
| 148 | +two_pi_f = 2 * np.pi * np.array([0.0, 2.0, 0.7]) |
| 149 | +two_pi_f_amp = np.multiply(two_pi_f, amp) |
132 | 150 | two_pi_f_squared_amp = np.multiply(two_pi_f, two_pi_f_amp) |
133 | 151 |
|
134 | 152 | for i in range(0, N_SIMULATION): |
135 | 153 | time_start = time.time() |
136 | | - |
137 | | - sampleCom.pos(offset + np.multiply(amp, np.sin(two_pi_f*t))) |
138 | | - sampleCom.vel(np.multiply(two_pi_f_amp, np.cos(two_pi_f*t))) |
139 | | - sampleCom.acc(np.multiply(two_pi_f_squared_amp, -np.sin(two_pi_f*t))) |
140 | | - |
| 154 | + |
| 155 | + sampleCom.pos(offset + np.multiply(amp, np.sin(two_pi_f * t))) |
| 156 | + sampleCom.vel(np.multiply(two_pi_f_amp, np.cos(two_pi_f * t))) |
| 157 | + sampleCom.acc(np.multiply(two_pi_f_squared_amp, -np.sin(two_pi_f * t))) |
| 158 | + |
141 | 159 | comTask.setReference(sampleCom) |
142 | 160 | samplePosture = trajPosture.computeNext() |
143 | 161 | postureTask.setReference(samplePosture) |
144 | 162 |
|
145 | 163 | HQPData = invdyn.computeProblemData(t, q, v) |
146 | | - if i == 0: HQPData.print_all() |
| 164 | + if i == 0: |
| 165 | + HQPData.print_all() |
147 | 166 |
|
148 | 167 | sol = solver.solve(HQPData) |
149 | | - if(sol.status!=0): |
150 | | - print("[%d] QP problem could not be solved! Error code:"%(i), sol.status) |
| 168 | + if sol.status != 0: |
| 169 | + print("[%d] QP problem could not be solved! Error code:" % (i), sol.status) |
151 | 170 | break |
152 | | - |
| 171 | + |
153 | 172 | tau = invdyn.getActuatorForces(sol) |
154 | 173 | dv = invdyn.getAccelerations(sol) |
155 | | - |
156 | | - com_pos[:,i] = robot.com(invdyn.data()) |
157 | | - com_vel[:,i] = robot.com_vel(invdyn.data()) |
158 | | - com_acc[:,i] = comTask.getAcceleration(dv) |
159 | | - com_pos_ref[:,i] = sampleCom.pos() |
160 | | - com_vel_ref[:,i] = sampleCom.vel() |
161 | | - com_acc_ref[:,i] = sampleCom.acc() |
162 | | - com_acc_des[:,i] = comTask.getDesiredAcceleration |
163 | | - |
164 | | - if i%PRINT_N == 0: |
165 | | - print("Time %.3f"%(t)) |
166 | | - print("\tNormal forces: ", end=' ') |
| 174 | + |
| 175 | + com_pos[:, i] = robot.com(invdyn.data()) |
| 176 | + com_vel[:, i] = robot.com_vel(invdyn.data()) |
| 177 | + com_acc[:, i] = comTask.getAcceleration(dv) |
| 178 | + com_pos_ref[:, i] = sampleCom.pos() |
| 179 | + com_vel_ref[:, i] = sampleCom.vel() |
| 180 | + com_acc_ref[:, i] = sampleCom.acc() |
| 181 | + com_acc_des[:, i] = comTask.getDesiredAcceleration |
| 182 | + |
| 183 | + if i % PRINT_N == 0: |
| 184 | + print("Time %.3f" % (t)) |
| 185 | + print("\tNormal forces: ", end=" ") |
167 | 186 | for contact in contacts: |
168 | 187 | if invdyn.checkContact(contact.name, sol): |
169 | 188 | f = invdyn.getContactForce(contact.name, sol) |
170 | | - print("%4.1f"%(contact.getNormalForce(f)), end=' ') |
| 189 | + print("%4.1f" % (contact.getNormalForce(f)), end=" ") |
171 | 190 |
|
172 | | - print("\n\ttracking err %s: %.3f"%(comTask.name.ljust(20,'.'), norm(comTask.position_error, 2))) |
173 | | - print("\t||v||: %.3f\t ||dv||: %.3f"%(norm(v, 2), norm(dv))) |
| 191 | + print( |
| 192 | + "\n\ttracking err %s: %.3f" |
| 193 | + % (comTask.name.ljust(20, "."), norm(comTask.position_error, 2)) |
| 194 | + ) |
| 195 | + print("\t||v||: %.3f\t ||dv||: %.3f" % (norm(v, 2), norm(dv))) |
174 | 196 |
|
175 | | - v_mean = v + 0.5*dt*dv |
176 | | - v += dt*dv |
177 | | - q = pin.integrate(robot.model(), q, dt*v_mean) |
| 197 | + v_mean = v + 0.5 * dt * dv |
| 198 | + v += dt * dv |
| 199 | + q = pin.integrate(robot.model(), q, dt * v_mean) |
178 | 200 | t += dt |
179 | | - |
180 | | - if i%DISPLAY_N == 0: robot_display.display(q) |
| 201 | + |
| 202 | + if i % DISPLAY_N == 0: |
| 203 | + robot_display.display(q) |
181 | 204 |
|
182 | 205 | time_spent = time.time() - time_start |
183 | | - if(time_spent < dt): time.sleep(dt-time_spent) |
| 206 | + if time_spent < dt: |
| 207 | + time.sleep(dt - time_spent) |
184 | 208 |
|
185 | 209 | # PLOT STUFF |
186 | | -time = np.arange(0.0, N_SIMULATION*dt, dt) |
| 210 | +time = np.arange(0.0, N_SIMULATION * dt, dt) |
187 | 211 |
|
188 | | -(f, ax) = plut.create_empty_figure(3,1) |
| 212 | +(f, ax) = plut.create_empty_figure(3, 1) |
189 | 213 | for i in range(3): |
190 | | - ax[i].plot(time, com_pos[i,:], label='CoM '+str(i)) |
191 | | - ax[i].plot(time, com_pos_ref[i,:], 'r:', label='CoM Ref '+str(i)) |
192 | | - ax[i].set_xlabel('Time [s]') |
193 | | - ax[i].set_ylabel('CoM [m]') |
| 214 | + ax[i].plot(time, com_pos[i, :], label="CoM " + str(i)) |
| 215 | + ax[i].plot(time, com_pos_ref[i, :], "r:", label="CoM Ref " + str(i)) |
| 216 | + ax[i].set_xlabel("Time [s]") |
| 217 | + ax[i].set_ylabel("CoM [m]") |
194 | 218 | leg = ax[i].legend() |
195 | 219 | leg.get_frame().set_alpha(0.5) |
196 | 220 |
|
197 | | -(f, ax) = plut.create_empty_figure(3,1) |
| 221 | +(f, ax) = plut.create_empty_figure(3, 1) |
198 | 222 | for i in range(3): |
199 | | - ax[i].plot(time, com_vel[i,:], label='CoM Vel '+str(i)) |
200 | | - ax[i].plot(time, com_vel_ref[i,:], 'r:', label='CoM Vel Ref '+str(i)) |
201 | | - ax[i].set_xlabel('Time [s]') |
202 | | - ax[i].set_ylabel('CoM Vel [m/s]') |
| 223 | + ax[i].plot(time, com_vel[i, :], label="CoM Vel " + str(i)) |
| 224 | + ax[i].plot(time, com_vel_ref[i, :], "r:", label="CoM Vel Ref " + str(i)) |
| 225 | + ax[i].set_xlabel("Time [s]") |
| 226 | + ax[i].set_ylabel("CoM Vel [m/s]") |
203 | 227 | leg = ax[i].legend() |
204 | 228 | leg.get_frame().set_alpha(0.5) |
205 | | - |
206 | | -(f, ax) = plut.create_empty_figure(3,1) |
| 229 | + |
| 230 | +(f, ax) = plut.create_empty_figure(3, 1) |
207 | 231 | for i in range(3): |
208 | | - ax[i].plot(time, com_acc[i,:], label='CoM Acc '+str(i)) |
209 | | - ax[i].plot(time, com_acc_ref[i,:], 'r:', label='CoM Acc Ref '+str(i)) |
210 | | - ax[i].plot(time, com_acc_des[i,:], 'g--', label='CoM Acc Des '+str(i)) |
211 | | - ax[i].set_xlabel('Time [s]') |
212 | | - ax[i].set_ylabel('CoM Acc [m/s^2]') |
| 232 | + ax[i].plot(time, com_acc[i, :], label="CoM Acc " + str(i)) |
| 233 | + ax[i].plot(time, com_acc_ref[i, :], "r:", label="CoM Acc Ref " + str(i)) |
| 234 | + ax[i].plot(time, com_acc_des[i, :], "g--", label="CoM Acc Des " + str(i)) |
| 235 | + ax[i].set_xlabel("Time [s]") |
| 236 | + ax[i].set_ylabel("CoM Acc [m/s^2]") |
213 | 237 | leg = ax[i].legend() |
214 | 238 | leg.get_frame().set_alpha(0.5) |
215 | | - |
| 239 | + |
216 | 240 | plt.show() |
0 commit comments