Parameter Identification of quarter car model

Hey guys,

I got the following system of a quarter car model:

The goal is to identify the unknown parameters of m1, m2, K, Kt and C based on differnt inputs q and outputs z1 and z2.

For this I simulated the model in Simulink with for example a sine wave as input q and set variables of m1, m2, K, Kt and C. I then recorded the position, motion and acceleration of z1 and z2.

Now lets say I only know the input q, output z1,z2 as well as their derivatives and I also got the equations of motion. I want to identify the unknown parameters m1, m2, K, Kt and C and see how close the identification gets compared to the values i defined in the simulink model.

I'd be glad if someone could give me hint as to how achieve this. Thanks for that.

[John D'Errico]

This is just a standard nonlinear least squares problem. You use an optimization tool, perhaps fminsearch, or one of the tools from the optimization or global optimization toolbox (or even the stats toolbox) to vary the unknown variables. Your objective function will need to run the simulation, taking the results, then return the error to the solver. Depending on the solver used, the objective function will be slightly different of course.

Answer 1

[Per Ljung]

This is a great homework problem. QCM is actually a linear simulation, and solving for the parameters is also linear so LS can be used.

The QCM has 2 masses and its behavior can be described using 4 generalized coordinates, for example z1, z2, z1d, z2d.

EOM

Balance the forces on each mass to derive the equations of motion as

% m1*z1dd = z1*(-k1) + z2*(k1) + z1d*(-c1) + z2d*(c1)
% m2*z2dd = z1*(k1)+z2*(-k1-k2)+z1d*(c1)+z2d*(-c1-c2)+q*(k2)

The above shows that you can compute the accelerations (z1dd,z2dd) from the states (z1,z2,z1d,z2d) and input (q).

Simulation

The QCM model is linear so it can be written using the matrices A,B,C,D and simulated. Here the C & D matrices include all the states (x) and accelerations in the output (y).

% QCM parameter values
m1=40; % masses [kg]
m2=5;
k2=150000; % spring rates [N/m]
k1=10000;
c1=550; % dampers [N*s/m]
c2=0;

V=5; % velocity [m/s]
T=0.5; % [sec]
dt=1e-3;

% road elevation profile q such as speedbump, curb, stochastic ...
% speedbump
BL=0.2; % bump length [m]
ds=V*dt; % 1cm
BS=[0:ds:BL] % bump distance [m]
BH=0.1; % bump height [m]
bump=BH*sin(pi*BS/BL);
% prepend & append flat sections of length L
L=1;
s=[0:ds:L, L+BS(2:end-1), L+BL+[0:ds:L]]' % want column
q=[0*[0:ds:L], bump(2:end-1), 0*[0:ds:L]]' % want column
t=s/V;

% state-space
A = [0,0,1,0; 
    0,0,0,1; 
    -k1/m1, k1/m1, -c1/m1, c1/m1;
    k1/m2, -(k1+k2)/m2, c1/m2, -(c1+c2)/m2 ];
B=[0; 0; 0; k2/m2];
C=[eye(4); 
    -k1/m1, k1/m1, -c1/m1, c1/m1; % z1dd
    k1/m2, -(k1+k2)/m2, c1/m2, -(c1+c2)/m2 ]; % z2dd
D=[0; 0; 0; k2/m2; 0; k2/m2];

% simulate
qcm=ss(A,B,C,D);
[y,tout,x]=lsim(qcm,q,t);

% states
z1=x(:,1);
z2=x(:,2);
z1d=x(:,3);
z2d=x(:,4);
% measurements
z1dd=y(:,5);
z2dd=y(:,6);

% plot results
figure
subplot(3,3,1)
plot(s,q,'k-'); title('q'); xlabel('dist [m]')
subplot(3,3,4)
plot(s,z1,'b-'); title('z1'); xlabel('dist [m]')
subplot(3,3,5)
plot(s,z1d,'b-'); title('z1d'); xlabel('dist [m]')
subplot(3,3,6)
plot(s,z1dd,'b-'); title('z1dd'); xlabel('dist [m]')
subplot(3,3,7)
plot(s,z2,'b-'); title('z2'); xlabel('dist [m]')
subplot(3,3,8)
plot(s,z2d,'b-'); title('z2d'); xlabel('dist [m]')
subplot(3,3,9)
plot(s,z2dd,'b-'); title('z2dd'); xlabel('dist [m]')
sgtitle("qcm CT simulation")

SysId

From the QCM simulation you now have the inputs (q) and outputs (z1,z1d,z1dd,z2,z2d,z2dd). Can you estimate m1,m2,c1,k1,k2 from the inputs and outputs? YES since it is a linear model, but not uniquely so you need some knowns e.g. z1dd,z2dd,m1,m2 which are easy to measure. Refactor the EOM to lhs=rhs*X where lhs are the knowns, rhs are the unknown factors, and X are the unknown parameters [k2; k1; c1]. Solve for X using over-determined least squares. This results in perfect parameter identification (assuming noiseless measurements).

% refactor eom with UNKNOWN parameters k2,k1,c1
% m1*z1dd = 0*(k2) -(z1-z2)*(k1) - (z1d-z2d)*(c1) 
% m2*z2dd = -(z2-q)*(k2) + (z1-z2)*(k1) + (z1d-z2d)*(c1)  
lhs = [m1*z1dd; m2*z2dd] 
rhs = [zeros(size(z1)), -(z1-z2), -(z1d-z2d); -(z2-q), (z1-z2), (z1d-z2d)]
% NOTE lhs and rhs must be column vectors

% solve for UNKNOWN parameters
X = rhs \ lhs % [k2; k1; c1] 
% 1.0e+05 * [1.5000;    0.1000;    0.0055]
% expected [k2;k1;c1]=[150000; 10000; 550]

% optional -- stdX is goodness indicator of each parameter estimate
[X,stdX,mse]=lscov(rhs,lhs)

% note that you can use less data (e.g. if you measured z2dd but not z1dd) 
% by just using the 2nd equation of motion
% m2*z2dd = - (z2-q)*k2 + (z1-z2)*k1 + (z1d-z2d)*c1  
lhs = [m2*z2dd] 
rhs = [-(z2-q), (z1-z2), (z1d-z2d)]
X = rhs \ lhs % [k2; k1; c1] 
% 1.0e+05 * [1.5000;    0.1000;    0.0055]
% expected [k2;k1;c1]=[150000; 10000; 550]