How to estimate state-space (ss) model of a mass-spring-damper system?

Dear all;

I have been trying to estimate the ss model of a mass-spring-damper system? it is 1-DOF very simple system and I assume that I already have the system parameters (mass, damping, stiffness, ...) and have simulated the response in time and frequency domains. Then I generated the model objects (frd, iddata). Interestingly Transfer-Function (TF) model is being estimated very good but state-space model always has isses and the system identification toolbox seems to be wrong.

I have both used the 'ident' (sys in App) and also the 'tfest' and 'ssest' built-in functions in the editor.

Any hints on improving the ss model estimation? (BTW, I have used different types of excitations signals: randn, sine, rbs, prbs, rgs, customized impulse, customized step, ... - code snippet is attaced)

% Plant Properties 
m1 = 1;
k1 = 1e6; 
ccr1 = 2*sqrt(m1*k1);
zeta1 = 0.1; 
c1 = ccr1 * zeta1; 
wn1 = sqrt(k1/m1);
wndam = wn1*sqrt(1-zeta1^2);
x10 = 0.25;
x10 = 0; % no IC meaning no static offset ... 
v10 = 0;
IC1 = [x10; v10];
% model the plant 
tfm1 = tf([0, 0, 1], [m1, c1, k1])
Am1 = [0, 1; 0-k1/m1, -c1/m1];
Bm1 = [0; 1/m1];
Cm1 = [1, 0];
Dm1 = 0;
ssm1 = ss(Am1, Bm1, Cm1, Dm1)
% Excitation / Driving Force / Input 
t1 = linspace(0, 2, 1000); % 1000 is L (total number of samples) 
u1 = 1 * sin(2*pi*t1); 
figure(1)
plot(t1, u1)
grid on  
% sys idn 
u2 = randn(size(t1)); 
[y1, tout1] = lsim(tfm1, u1, t1, IC1); % lsim simulates the sys response 
[y2, tout2] = lsim(tfm1, u2, t1, IC1); 
u3 = ones(size(t1));
[y3, tout3] = lsim(tfm1, u3, t1, IC1);
u4 = zeros(size(t1));
u4(t1 <= 0.01) = 1;  
[y4, tout4] = lsim(tfm1, u4, t1, IC1); %impulse(tfm1);
% [y5, tout5] = ramp(tfm1);
[H1, wout1] = freqresp(ssm1); 
figure(2)
subplot(2, 1, 1)
plot(t1', u1');
grid on  
subplot(2, 1, 2)
plot(tout1, y1)
grid on 
figure(3)
subplot(2, 1, 1)
plot(t1', u2');
grid on  
subplot(2, 1, 2)
plot(tout1, y2)
grid on 
figure(4)
step(tfm1)
grid on
ts = 0.002; % 
% Defining Model Objects (specialty data container)
mo1 = frd(H1, wout1); % Model Object 1 

mo21 = iddata(y1, u1', ts); % Model Object 2
mo22 = iddata(y2, u2', ts); % 
mo23 = iddata(y3, u3', ts); 
mo24 = iddata(y4, u4', ts); % 
% freq. resp.
% Perform System Identification 
nz = 0; % number of zeros 
np = 2; % number of poles 
nx = 2; % number of states 
% [Estsys1, ic1] = ssest(mo1, nx)
[Esttfmo1, ictfmo1] = tfest(mo1, np, nz) % ic1 is the estimated IC 
[partfmo1, stdtfmo1] = getpvec(Esttfmo1)
Esttfmo1cov = getcov(Esttfmo1)
[Estssmo1, icssmo1] = ssest(mo1, nx) % ic1 is the estimated IC 
[parssmo1, stdssmo1] = getpvec(Estssmo1)
Estssmo1cov = getcov(Estssmo1)
% use iddata mo2... - tfest 
[Esttfmo21, ictfmo21] = tfest(mo21, np, nz)
[partfmo21, stdtfmo21] = getpvec(Esttfmo21)
[Esttfmo22, ictfmo22] = tfest(mo22, np, nz)
[partfmo22, stdtfmo22] = getpvec(Esttfmo22)
[Esttfmo23, ictfmo23] = tfest(mo23, np, nz)
[partfmo23, stdtfmo23] = getpvec(Esttfmo23)
[Esttfmo24, ictfmo24] = tfest(mo24, np, nz)
[partfmo24, stdtfmo24] = getpvec(Esttfmo24)
% use iddata mo2... - ssest 
[Estssmo21, icssmo21] = ssest(mo21, nx)
[parssmo21, stdssmo21] = getpvec(Estssmo21)
[Estssmo22, icssmo22] = ssest(mo22, nx)
[parssmo22, stdssmo22] = getpvec(Estssmo22)
[Estssmo23, icssmo23] = ssest(mo23, nx)
[parssmo23, stdssmo23] = getpvec(Estssmo23)
[Estssmo24, icssmo24] = ssest(mo24, nx)
[parssmo24, stdssmo24] = getpvec(Estssmo24)
%% Converting  aContinous-time transfer function into Discrete-time tf via: c2d
tfm1d = c2d(tfm1, 0.01); % 0.01 here is the sampling time 
disp(tfm1)
disp(tfm1d)
report1 = Esttfmo1.Report
% report1_ics = Estsys1.Report.InitialState
figure(6)
stepplot(tfm1, '-b')
hold on 
stepplot(Esttfmo1, '*r')
hold on 
% stepplot(Estsys2, '-.k')
% hold on 
% stepplot(Estsys3, '*y')
% grid on 
legend({'Simulation Response ro Step Excitations', ...
    'Estimated Model Response to Step Excitations (tfest(idfrd))', ...
    'Estimated Model Response to Step Excitations (tfest(iddata))', ...
       'Estimated Model Response to Step Excitations (ssest(iddata))'}, ...
    'FontSize', 14)

figure(7)
subplot(2, 1, 1)
plot(wout1, abs(squeeze(H1(1, 1, :))))
grid on 
subplot(2, 1, 2)
plot(wout1, abs(squeeze(H1(1, 1, :))))
grid on 

figure(8)
subplot(2, 1, 1)
loglog(wout1, abs(squeeze(H1(1, 1, :))))
grid on 
subplot(2, 1, 2)
loglog(wout1, 20*abs(squeeze(H1(1, 1, :))))
grid on 

noise1 = randn(2*size(wout1));
noise11 = fft(noise1); 
nfft1 = abs(noise11 ./ (size(noise1)));
nfft2 = nfft1(1 : (size(noise1)/2)+1);
nfft2(2 : end-1) = 2 * nfft2(2 : end-1);
nfft3 = nfft2(1 : end-1); 
nfft3 = 1e-9 .* nfft3;

% H1sonly = squeeze(H1); % to remove the ingeleton!!!
H1only = reshape(H1, [], 1);  
H1noise2 = randn(size(H1only)); 
noise2 = 1e-6 .* H1noise2;
noise3 = abs(noise2); 
size(H1only) 
H1n1 = H1only + nfft3'; 
H1n1 = H1only + noise3; 
size(H1n1)
figure(9)
loglog(wout1, abs(squeeze(H1(1, 1, :))))
grid on 
hold on 
loglog(wout1, H1n1)
figure(10)
compare(mo1, Esttfmo1)
figure(11)
compare(mo1, Estssmo1)
figure(12)
compare(mo1, Esttfmo24)
figure(13)
compare(mo1, Estssmo23)
figure(14)
compare(mo1, Estssmo24)

figure(16)
plot(tout2, y2)

% so far: impulse, step, randn did NOT work in terms of sys idn. 
NN = 1000;%200; % number of samples in the excitation signal 
u11 = idinput(NN, 'rbs'); % random binary signal - values either -1 or +1 
% mo31 = iddata([], u11)
mo31 = iddata([], u11, 1)
figure(31)
plot(u11)
grid on 
title('rbs: Random Binary Signal (-1 or +1)', 'FontSize', 22)
u12 = idinput(NN, 'prbs'); % Pseudo-random binary signal - values either -1 or +1 
mo32 = iddata([], u12, 1)
figure(32)
plot(u12)
grid on 
title('prbs: PseudoRandom Binary Signal', 'FontSize', 22)
u13 = idinput(NN, 'rgs'); % Random Gaussian signal - values not deterministic 
mo33 = iddata([], u13, 1)
figure(33)
plot(u13)
grid on 
title('rgs: Random Gaussian Signal', 'FontSize', 22)
u14 = idinput(NN, 'sine'); % Multi-Sine  
mo34 = iddata([], u14, 1)
figure(34)
plot(u14)
grid on 
title('sine: Multi-Sine Signal', 'FontSize', 22)
% 
[y31, tout31] = lsim(ssm1, u11, t1, IC1)
iddata11 = iddata(y31, u11, 0.001);
[y32, tout32] = lsim(ssm1, u12, t1, IC1)
iddata12 = iddata(y32, u12, 0.001);
[y33, tout33] = lsim(ssm1, u13, t1, IC1)
iddata13 = iddata(y33, u13, 0.001);
[y34, tout34] = lsim(ssm1, u14, t1, IC1)
iddata14 = iddata(y34, u14, 0.001);
% sys idn 
etf11 = tfest(iddata11, np, nz);
etf12 = tfest(iddata12, np, nz);
etf13 = tfest(iddata13, np, nz);
etf14 = tfest(iddata14, np, nz);

ess11 = ssest(iddata11, nx);
ess12 = ssest(iddata12, nx);
ess13 = ssest(iddata13, nx);
ess14 = ssest(iddata14, nx);

Thanks

[Tianyu]

You have a long script, can you please specify what issue did you encounter? In which line? What is your expected behavior?

The best practice is to provide a shorter code only focusing on the issue you mentioned.

[Mathieu NOE]

hello @Alireza

I don't use the system identification toolbox... maybe because I am still using mostly tools from the past that are still very valuable .

have you tried the subspace identification routines created by Peter Van Overschee ?

attached is a zip of his tollbox (main routine is subid.m)

hope it helps !

[Mathieu NOE]

his work