Sebastian Uribe / Mbed OS pan_flipping

Dependencies:   MatrixMath Matrix ExperimentServer QEI_pmw MotorShield

Revision:
34:9a24d0f718ac
Parent:
32:c60a5d33cd79
Child:
35:88dbfefc1bbb
--- a/main.cpp	Mon Nov 16 05:12:37 2020 +0000
+++ b/main.cpp	Mon Nov 16 08:23:00 2020 +0000
@@ -6,11 +6,9 @@
 #include "BezierCurve.h"
 #include "MotorShield.h" 
 #include "HardwareSetup.h"
-#include "Matrix.h"
-#include "MatrixMath.h"
 
 #define BEZIER_ORDER_FOOT    7
-#define NUM_INPUTS (12 + 2*(BEZIER_ORDER_FOOT+1))
+#define NUM_INPUTS (14 + 2*(BEZIER_ORDER_FOOT+1))
 #define NUM_OUTPUTS 19
 
 #define PULSE_TO_RAD (2.0f*3.14159f / 1200.0f)
@@ -28,20 +26,15 @@
 MotorShield motorShield(12000); //initialize the motor shield with a period of 12000 ticks or ~20kHZ
 Ticker currentLoop;
 
-Matrix MassMatrix(2,2);
-Matrix Jacobian(2,2);
-Matrix JacobianT(2,2);
-Matrix InverseMassMatrix(2,2);
-Matrix temp_product(2,2);
-Matrix Lambda(2,2);
-
 // Variables for q1
 float current1;
 float current_des1 = 0;
 float prev_current_des1 = 0;
 float current_int1 = 0;
 float angle1;
+float angle_des1;
 float velocity1;
+float velocity_des1;
 float duty_cycle1;
 float angle1_init;
 
@@ -51,32 +44,26 @@
 float prev_current_des2 = 0;
 float current_int2 = 0;
 float angle2;
+float angle_des2;
 float velocity2;
+float velocity_des2;
 float duty_cycle2;
 float angle2_init;
 
+// Fixed kinematic parameters
+const float l_OA=.011; 
+const float l_OB=.042; 
+const float l_AC=.096; 
+const float l_DE=.091;
 
+float l1; //=0.084125;
+float l2; //=0.084125;
 
 // Timing parameters
 float current_control_period_us = 200.0f;     // 5kHz current control loop
 float impedance_control_period_us = 1000.0f;  // 1kHz impedance control loop
 float start_period, traj_period, end_period;
 
-// Hardware kinematic parameters -- NEED PAN PARAMETERS
-const float l_c1; //upper arm center of mass
-const float l_B; //upper arm length
-const float r_c2; //lower arm center of mass
-const float l_C; //lower arm length
-
-const float m1; //mass of upper arm
-const float m2; //mass of lower arm
-
-const float I1; //upper arm interia
-const float I2;  //lower arm inertia
-
-const float N; //gear ratio
-const float Ir; //motor inertia
-
 // Control parameters
 float current_Kp = 4.0f;         
 float current_Ki = 0.4f;           
@@ -103,7 +90,7 @@
     //use the motor shield as follows:
     //motorShield.motorAWrite(DUTY CYCLE, DIRECTION), DIRECTION = 0 is forward, DIRECTION =1 is backwards.
         
-    current1 = -(((float(motorShield.readCurrentA())/65536.0f)*30.0f)-15.0f);           // measure current
+    current1 = -(((float(motorShield.readCurrentA())/65536.0f)*30.0f)-15.0f);           // measure current DO WE NEED TO ADJUST GEAR RATIO HERE?
     velocity1 = encoderA.getVelocity() * PULSE_TO_RAD;                                  // measure velocity        
     float err_c1 = current_des1 - current1;                                             // current errror
     current_int1 += err_c1;                                                             // integrate error
@@ -123,7 +110,7 @@
     }             
     prev_current_des1 = current_des1; 
     
-    current2     = -(((float(motorShield.readCurrentB())/65536.0f)*30.0f)-15.0f);       // measure current
+    current2     = -(((float(motorShield.readCurrentB())/65536.0f)*30.0f)-15.0f);       // measure current 
     velocity2 = encoderB.getVelocity() * PULSE_TO_RAD;                                  // measure velocity  
     float err_c2 = current_des2 - current2;                                             // current error
     current_int2 += err_c2;                                                             // integrate error
@@ -142,10 +129,12 @@
         motorShield.motorBWrite(absDuty2, 0);
     }             
     prev_current_des2 = current_des2; 
+    
 }
 
 int main (void)
 {
+    
     // Object for 7th order Cartesian foot trajectory
     BezierCurve rDesFoot_bez(2,BEZIER_ORDER_FOOT);
     
@@ -158,37 +147,28 @@
     pc.printf("%f",input_params[0]);
     
     while(1) {
+        
         // If there are new inputs, this code will run
         if (server.getParams(input_params,NUM_INPUTS)) {
+            
+                        
             // Get inputs from MATLAB          
-            start_period = input_params[0]; // First buffer time, before trajectory
-            traj_period = input_params[1]; // Trajectory time/length
-            end_period = input_params[2]; // Second buffer time, after trajectory
-            
-            l_c1 = input_params[3]; //upper arm center of mass
-            l_B = input_params[4]; //upper arm length
-            l_c2 = input_params[5]; //lower arm center of mass
-            l_C = input_params[6]; //lower arm length
-            
-            m1 = input_params[7]; //mass of upper arm
-            m2 = input_params[8]; //mass of lower arm
-            
-            I1 = input_params[9]; //upper arm interia
-            I2 = input_params[10];  //lower arm inertia
-            
-            N = input_params[11]; //gear ratio
-            Ir = input_params[12]; //motor inertia
+            start_period                = input_params[0];    // First buffer time, before trajectory
+            traj_period                 = input_params[1];    // Trajectory time/length
+            end_period                  = input_params[2];    // Second buffer time, after trajectory
     
-            angle1_init = input_params[13]; // Initial angle for q1 (rad)
-            angle2_init = input_params[14]; // Initial angle for q2 (rad)
+            angle1_init                 = input_params[3];    // Initial angle for q1 (rad)
+            angle2_init                 = input_params[4];    // Initial angle for q2 (rad)
 
-            K_xx = input_params[15]; // Foot stiffness N/m
-            K_yy = input_params[16]; // Foot stiffness N/m
-            K_xy = input_params[17]; // Foot stiffness N/m
-            D_xx = input_params[18]; // Foot damping N/(m/s)
-            D_yy = input_params[19]; // Foot damping N/(m/s)
-            D_xy = input_params[20]; // Foot damping N/(m/s)
-            duty_max = input_params[21]; // Maximum duty factor
+            K_xx                        = input_params[5];    // Foot stiffness N/m
+            K_yy                        = input_params[6];    // Foot stiffness N/m
+            K_xy                        = input_params[7];    // Foot stiffness N/m
+            D_xx                        = input_params[8];    // Foot damping N/(m/s)
+            D_yy                        = input_params[9];    // Foot damping N/(m/s)
+            D_xy                        = input_params[10];   // Foot damping N/(m/s)
+            duty_max                    = input_params[11];   // Maximum duty factor
+            l1 = input_params[12];
+            l2 = input_params[13];
           
             // Get foot trajectory points
             float foot_pts[2*(BEZIER_ORDER_FOOT+1)];
@@ -227,27 +207,35 @@
                 const float dth2= velocity2;
  
                 // Calculate the Jacobian
-                float Jx_th1 = l_AC*cos(th1+th2)+l_DE*cos(th1)+l_OB*cos(th1);
-                float Jx_th2 = l_AC*cos(th1+th2);
-                float Jy_th1 = l_AC*sin(th1+th2)+l_DE*sin(th1)+l_OB*sin(th1);
-                float Jy_th2 = l_AC*sin(th1+th2);
-                                
+                float Jx_th1 = l1*cos(th1)+l2*cos(th1+th2);
+                float Jx_th2 = l2*cos(th1+th2);
+                float Jy_th1 = l1*sin(th1)+l2*sin(th1+th2);
+                float Jy_th2 = l2*sin(th1+th2);
+                
+
                 // Calculate the forward kinematics (position and velocity)
-                float xFoot = l_DE*sin(th1)+l_OB*sin(th1)+l_AC*sin(th1+th2);
-                float yFoot = -l_DE*cos(th1)-l_OB*cos(th1)-l_AC*cos(th1+th2);
-                float dxFoot = dth1*(l_AC*cos(th1+th2)+l_DE*cos(th1)+l_OB*cos(th1))+dth2*l_AC*cos(th1+th2);
-                float dyFoot = dth1*(l_AC*sin(th1+th2)+l_DE*sin(th1)+l_OB*sin(th1))+dth2*l_AC*sin(th1+th2); 
- 
-
+//                float xFoot = l_DE*sin(th1)+l_OB*sin(th1)+l_AC*sin(th1+th2);
+//                float yFoot = -l_DE*cos(th1)-l_OB*cos(th1)-l_AC*cos(th1+th2);
+//                float dxFoot = dth1*(l_AC*cos(th1+th2)+l_DE*cos(th1)+l_OB*cos(th1))+dth2*l_AC*cos(th1+th2);
+//                float dyFoot = dth1*(l_AC*sin(th1+th2)+l_DE*sin(th1)+l_OB*sin(th1))+dth2*l_AC*sin(th1+th2);  
+                float xFoot = l2*sin(th1+th2) + l1*sin(th1);
+                float yFoot = -l2*cos(th1+th2) - l1*cos(th1);
+                float dxFoot = l1*cos(th1)*dth1 + l2*cos(th1+th2)*(dth1+dth2);
+                float dyFoot = l1*sin(th1)*dth1 + l2*sin(th1+th2)*(dth1+dth2);
+                
                 // Set gains based on buffer and traj times, then calculate desired x,y from Bezier trajectory at current time if necessary
                 float teff  = 0;
                 float vMult = 0;
                 if( t < start_period) {
-                    if (K_xx > 0 || K_yy > 0) {
-                        K_xx = 100; 
-                        K_yy = 100; 
-                        D_xx = 5;  
-                        D_yy = 5;  
+                        if (K_xx > 0 || K_yy > 0) {
+//                        K_xx = 1; // for joint space control, set this to 1
+//                        K_yy = 1; // for joint space control, set this to 1
+//                        D_xx = .1;  // for joint space control, set this to 0.1
+//                        D_yy = .1;  // for joint space control, set this to 0.1
+                        K_xx = 50; // for joint space control, set this to 1
+                        K_yy = 50; // for joint space control, set this to 1
+                        D_xx = 2;  // for joint space control, set this to 0.1
+                        D_yy = 2;  // for joint space control, set this to 0.1
                         K_xy = 0;
                         D_xy = 0;
                     }
@@ -255,12 +243,12 @@
                 }
                 else if (t < start_period + traj_period)
                 {
-                    K_xx = input_params[15];  // Foot stiffness N/m
-                    K_yy = input_params[16];  // Foot stiffness N/m
-                    K_xy = input_params[17];  // Foot stiffness N/m
-                    D_xx = input_params[18];  // Foot damping N/(m/s)
-                    D_yy = input_params[19];  // Foot damping N/(m/s)
-                    D_xy = input_params[20]; // Foot damping N/(m/s)
+                    K_xx = input_params[5];  // Foot stiffness N/m
+                    K_yy = input_params[6];  // Foot stiffness N/m
+                    K_xy = input_params[7];  // Foot stiffness N/m
+                    D_xx = input_params[8];  // Foot damping N/(m/s)
+                    D_yy = input_params[9];  // Foot damping N/(m/s)
+                    D_xy = input_params[10]; // Foot damping N/(m/s)
                     teff = (t-start_period);
                     vMult = 1;
                 }
@@ -270,7 +258,6 @@
                     vMult = 0;
                 }
                 
-                // Get desired foot positions and velocities
                 float rDesFoot[2] , vDesFoot[2];
                 rDesFoot_bez.evaluate(teff/traj_period,rDesFoot);
                 rDesFoot_bez.evaluateDerivative(teff/traj_period,vDesFoot);
@@ -292,57 +279,34 @@
                 float dth2_des = (1.0f/dd) * ( -Jy_th1*vDesFoot[0] + Jx_th1*vDesFoot[1] );
         
                 // Calculate error variables
-                float e_x = 0;
-                float e_y = 0;
-                float de_x = 0;
-                float de_y = 0;
+                float e_x = xFootd-xFoot;
+                float e_y = yFootd-yFoot;
+                float de_x = vDesFoot[0]-dxFoot;
+                float de_y = vDesFoot[1]-dyFoot;
         
                 // Calculate virtual force on foot
-                float xdelta = -xFootd - xFoot;
-                float ydelta = yFootd - yFoot;
-                float dydelta = vDesFoot[1] - dyFoot;
-                float dxdelta = vDesFoot[0] - dxFoot;
-                float fx = K_xx*xdelta+K_xy*ydelta+D_xx*dxdelta+D_xy*dydelta;
-                float fy = K_yy*ydelta+K_xy*xdelta+D_yy*dydelta+D_xy*dxdelta;
-                
-                float t1 = Jx_th1*fx + Jy_th1*fy;
-                float t2 = Jx_th2*fx + Jy_th2*fy;
-                
-                // Calculate mass matrix elements
-                float M11 = I1 + I2 + I3 + I4 + Ir + Ir*N*N + l_AC*l_AC*m4 + l_A_m3*l_A_m3*m3 + l_B_m2*l_B_m2*m2 + l_C_m4*l_C_m4*m4 + l_OA*l_OA*m3 + l_OB*l_OB*m2 + l_OA*l_OA*m4 + l_O_m1*l_O_m1*m1 + 2*l_C_m4*l_OA*m4 + 2*l_AC*l_C_m4*m4*cos(th2) + 2*l_AC*l_OA*m4*cos(th2) + 2*l_A_m3*l_OA*m3*cos(th2) + 2*l_B_m2*l_OB*m2*cos(th2); 
-                float M12 = I2 + I3 + l_AC*l_AC*m4 + l_A_m3*l_A_m3*m3 + l_B_m2*l_B_m2*m2 + Ir*N + l_AC*l_C_m4*m4*cos(th2) + l_AC*l_OA*m4*cos(th2) + l_A_m3*l_OA*m3*cos(th2) + l_B_m2*l_OB*m2*cos(th2); 
-                float M22 = Ir*N*N + m4*l_AC*l_AC + m3*l_A_m3*l_A_m3 + m2*l_B_m2*l_B_m2 + I2 + I3;
- 
-                // Populate mass matrix
-                MassMatrix.Clear();
-                MassMatrix << M11 << M12
-                           << M12 << M22;
+                float fx = K_xx*e_x + K_xy*e_y + D_xx*de_x + D_xy*de_y;
+                float fy = K_xy*e_x + K_yy*e_y + D_xy*de_x + D_yy*de_y;
                 
-                // Populate Jacobian matrix
-                Jacobian.Clear();
-                Jacobian << Jx_th1 << Jx_th2
-                         << Jy_th1 << Jy_th2;
-                
-                // Calculate Lambda matrix
-                JacobianT = MatrixMath::Transpose(Jacobian);
-                InverseMassMatrix = MatrixMath::Inv(MassMatrix);
-                temp_product = Jacobian*InverseMassMatrix*JacobianT;
-                Lambda = MatrixMath::Inv(temp_product); 
-                
-                // Pull elements of Lambda matrix
-                float L11 = Lambda.getNumber(1,1);
-                float L12 = Lambda.getNumber(1,2);
-                float L21 = Lambda.getNumber(2,1);
-                float L22 = Lambda.getNumber(2,2);      
-                
-                // Calculate desired motor torque
-                float t1_op = (Jx_th1*L11+Jy_th1*L21)*fx + (Jx_th1*L12+Jy_th1*L22)*fy;
-                float t2_op = (Jx_th2*L11+Jy_th2*L21)*fx + (Jx_th2*L12+Jy_th2*L22)*fy;      
+                float T1 = Jx_th1*fx+Jy_th1*fy;    
+                float T2 = Jx_th2*fx+Jy_th2*fy;    
+                    
                                 
                 // Set desired currents             
-                current_des1 = t1_op/k_t;          
-                current_des2 = t2_op/k_t;   
-
+        
+                // Joint impedance
+                // sub Kxx for K1, Dxx for D1, Kyy for K2, Dyy for D2
+                // Note: Be careful with signs now that you have non-zero desired angles!
+                // Your equations should be of the form i_d = K1*(q1_d - q1) + D1*(dq1_d - dq1)
+//                current_des1 = (K_xx*(th1_des-angle1) + D_xx*(dth1_des-velocity1))/k_t;          
+//                current_des2 = (K_yy*(th2_des-angle2) + D_yy*(dth2_des-velocity2))/k_t;                       
+                           
+                // Cartesian impedance  
+                // Note: As with the joint space laws, be careful with signs!              
+                current_des1 = T1/k_t;          
+                current_des2 = T2/k_t;
+                
+                
                 // Form output to send to MATLAB     
                 float output_data[NUM_OUTPUTS];
                 // current time
@@ -380,6 +344,9 @@
             currentLoop.detach();
             motorShield.motorAWrite(0, 0); //turn motor A off
             motorShield.motorBWrite(0, 0); //turn motor B off
+        
         } // end if
+        
     } // end while
+    
 } // end main
\ No newline at end of file