least square method with multible unkonwn

Hello Everyone,

I am trying to fit some equations to data. As I am quite new to Matlab and coding in general, I am having a hard time achieving this.

Maybe you can help me out with my problem.

I got data from three 3-axis-accelerometers and I am trying to get the angular velocity/acceleration.

The vehicle was not moving during testing, therfore it should be feasible to say that the global coordinate system is equal to the local coordinate system. If errors are very high at the end, I might have to change this.

The measured data looks like this:

%Time	Acceleration	Acceleration	Acceleration	Acceleration	Acceleration	Acceleration	Acceleration	Acceleration	Acceleration
%s	m/s^2	m/s^2	m/s^2	m/s^2	m/s^2	m/s^2	m/s^2	m/s^2	m/s^2
%	CylHead2x 	CylHead2y 	CylHead2z 	CylHead1x	CylHead1y	CylHead1z	MAGx	MAGy	MAGz
0	0	0	0	0	0	0	0	0	0
2,08333E-05	6,33128E-08	-4,59763E-08	2,73223E-08	-5,05531E-08	-2,49645E-08	-4,32718E-08	-4,84728E-08	-6,31047E-09	-3,13443E-08
4,16667E-05	5,5274E-07	-3,48914E-07	2,55787E-07	-5,98569E-09	-7,69206E-08	-1,61455E-07	-3,42487E-07	7,88509E-09	-3,88854E-07
0,0000625	2,44903E-06	-1,14245E-06	1,06341E-06	1,79757E-06	-2,49593E-07	-3,15895E-07	-1,34916E-06	1,37689E-07	-1,90947E-06
8,33333E-05	7,35567E-06	-2,01004E-06	2,81663E-06	9,60767E-06	5,01715E-07	5,32651E-07	-3,46298E-06	3,39389E-07	-5,71277E-06
0,000104167	1,71466E-05	-1,31856E-06	5,2826E-06	2,83488E-05	2,21753E-06	8,64681E-06	-6,71597E-06	3,42216E-07	-1,23976E-05
0,000125	3,42688E-05	2,97117E-06	6,73722E-06	5,97391E-05	2,41451E-06	3,5692E-05	-1,2059E-05	2,0028E-07	-2,10872E-05
%...

First I read the data which worked fine:

filename = 'a.xlsx';
sheet = 1;
xlRange = 'A4:J10000';
A = xlsread(filename,sheet,xlRange);

time    = readvars('a.xlsx','sheet','Sheet1','Range','A4:A50');
cyl1x   = readvars('a.xlsx','sheet','Sheet1','Range','E4:E50');
cyl1y   = readvars('a.xlsx','sheet','Sheet1','Range','F4:F50');
cyl1z   = readvars('a.xlsx','sheet','Sheet1','Range','G4:G50');
cyl2x   = readvars('a.xlsx','sheet','Sheet1','Range','B4:B50');
cyl2y   = readvars('a.xlsx','sheet','Sheet1','Range','C4:C50');
cyl2z   = readvars('a.xlsx','sheet','Sheet1','Range','D4:D50');
MAGx    = readvars('a.xlsx','sheet','Sheet1','Range','H4:H50');
MAGy    = readvars('a.xlsx','sheet','Sheet1','Range','I4:I50');
MAGz    = readvars('a.xlsx','sheet','Sheet1','Range','J4:J50');

%distance from accelerometers to CoG
r_z1 = [   27     268      39];   % Cyl 1  
r_z2 = [   58    -260      39];   % Cyl 2  
r_MAG= [ -279     -58      27];   % MAG 

Next I tried to get solutions for the following dependend equations in the least square sense:

Unknown: wp = angular acceleration[edited], w = angular velocity

cyl1x = cyl2x+wpy(r_z1(3)-r_z2(3))-wpz*(r_z1(2)-r_z2(2))+wx*(wy*(r_z1(2)-r_z2(2))+wz*(r_z1(3)-r_z2(3)))-(wy^2+wz^2)*(r_z1(1)-r_z2(1));
cyl1y = MAGy+wpz(r_z1(1)-r_MAG(1))-wpx*(r_z1(3)-r_MAG(3))+wy*(wx*(r_z1(1)-r_MAG(1))+wz*(r_z1(3)-r_MAG(3)))-(wx^2+wz^2)*(r_z1(2)-r_MAG(2));
cyl1z = cyl2z+wpx(r_z1(2)-r_z2(2))-wpy*(r_z1(1)-r_z2(1))+wz*(wx*(r_z1(1)-r_z2(1))+wy*(r_z1(2)-r_z2(2)))-(wx^2+wy^2)*(r_z1(3)-r_z2(3));

I tried to solve this with a time loop and linear least square method. Unfortunatly I am having a hard time getting anything usefull here. If you can give me a hint how to do this right I would be really greatful.

Best regards, Rolfe

[Torsten]

Don't you have to use somehow that angular velocity is the integral of angular acceleration with respect to time ?

And do I get it right that for each time instant, you have three equations for six unknowns ?

And why do you overwrite cyl1x every time when reading from the Excel file ?

cyl1x   = readvars('a.xlsx','sheet','Sheet1','Range','E4:E50');
cyl1x   = readvars('a.xlsx','sheet','Sheet1','Range','F4:F50');
cyl1x   = readvars('a.xlsx','sheet','Sheet1','Range','G4:G50');

[Rolfe]

Yes i do need the integral for the angular velocity but i dont know how to include that in the least square method. With this there are only 3 unknown and 3 equations right?

Sorry in my code I didn't overwrite it. Must have been a mistake when writing this question. This should just be Cyl1x, Cyl1y, Cyl1z

[Torsten]

With this there are only 3 unknown and 3 equations right?

Right. By integrating the acceleration data, you should get estimates for velocity. But if you have measured acceleration, why do you try to calculate it from your equations ? Or is the calculated acceleration different from the measured acceleration ?

[Rolfe]

I need the angular velocity and the angular acceleration of the body. These are my unknowns.

My measured data are only linear acceleration of different point (not CoG) of that body. To get the acceleration in CoG i need the angular velocity as well as the angular acceleration. With these equations I try to get the unknown.

[Torsten]

Ok, then you will have to solve the equations for at least two time instants together and express the angular velocity at time t as a function of the angular acceleration at times t and t-dt (or something similar - I'm not an expert in mechanics).

[Rolfe]

Yes, this is exactly where I am stuck. Not the mechanics behind it but the corresponding matlabcode to get to the right solution.

[Torsten]

It's not a MATLAB question.

It's a question on how usually wx(i), wy(i) and wz(i) are approximated as functions of w.x(i),w.y(i), w.z(i) and maybe w.x(i-1),w.y(i-1), w.z(i-1),w.x(i+1),w.y(i+1), w.z(i+1).

Answer 1

[William Rose]

@Rolfe, I will look at this later this week, if you do not have an answer by then.

Below is a part of an analysis I did related to soccer heading data. I have to think about this more. I will have time later this week.

We will assume that the center of mass is point A, and we will assume that the linear acceleration of the center of mass is zero: . In the equaiton below, we will use B=B1, B2, B3 for the 2 cylinders and the magneto respectively. We will try to find values for ω and that minimize the sum squared error between the estimated and measured linear accelearations at the three points.

I hope the solution does not involve a Kalman filter. I never mastered Kalman filters.

1. Angular velocity estimation based on adaptive simplified spherical simplex unscented Kalman filter in GFSINS
2. Kalman Filter estimation of angular acceleration https://iopscience.iop.org/article/10.1088/1757-899X/916/1/012072/pdf

----------------

All points on a rigid body experience the same angular velocity and angular acceleration at each time point. Let A and B be two points on a rigid body. The linear acceleration at point B is related to the movement at point A by

where

= linear accelerations at points A and B

ω = angular velocity (same everywhere, so it is not indexed by A or B)

α = dω/dt = angular acceleration (same everywhere, so it is not indexed by A or B)

= vector pointing to B from A.

Sources: “Chapter 15 Kinematics of Rigid Bodies”, retrieved 2016-01-27 from here, see equation near top of p.23. Whiteman, W, “Planar (2D) Rigid Body Kinematics I”, retrieved 2016-01-27 from here, see time 2:40. Shaprio, J, “Chapter 4 Rigid Body Motion”, retrieved 20160128 from here; see equation at bottom of p.12; this equation includes two extra terms which drop out, since, for a rigid body, the intra-body velocity db/dt and acceleration are 0.

[Rolfe]

Thank you so much for your answer William!

The three equations on top are derived from the equation you mentioned. I had 6 equations (ap1x, ap1y, ap1z, ap2x, ap2z, ap3y) and with these i eliminated the acceleration of the center of mass. So in the equation above aB and aA are both measured points (Cyl1, Cyl2, MAG).

"We will assume that the center of mass is point A, and we will assume that the linear acceleration of the center of mass is zero" I think I can't say that the linear acceleration of the center of mass is zero as this is part of the solution i would need. Or did I understand that wrong?

[William Rose]

@Rolfe,

[edit: correct typos]

OK, let us not assume the center of mass linear acceleration is zero. Let us try to solve for the CoM linear acceleration, using the data available.

You have

%distance from accelerometers to CoG
r_z1 = [   27     268      39];   % Cyl 1  
r_z2 = [   58    -260      39];   % Cyl 2  
r_MAG= [ -279     -58      27];   % MAG 

I assume the units above are mm, for an engine of reasonable size. The units for your linear acceleration data are m/s^2. Use the same distance units for both. This will not solve the problem, but will be important sometime.

Can you post a file of acceleration data? I see that the sampling rate is 48 kHz.

[Rolfe]

@William Rose Yes the units for the distance are in mm. I will change that right away.

"Can you post a file of acceleration data?"

sure thing. See attached file.

[William Rose]

@Rolfe, thank you. See email.