estimate(system)

function P = estimate(system)

% Time data
tsamp = 0.001;
t = (0:tsamp:6.0)';

% Create Input data (normally this come from some experiment)
U = 0.06;
fi =  0.1;
ff = 20.0;
tf =  6.0;
u = U*sin((fi + (ff-fi)/tf*t) .* t);

% Nominal model parameters
cv = 0.47;
fd = 4.7;
m  = 1/4.7;
ks = 1220;
var = 0.0005; % Artificial noise
exactP = [cv, fd];

fixParam = {'m', 'ks', 'var'}; % Fixed model parameter
fixParamNum = [m, ks, var];
varParam = {'cv', 'fd'};       % Unknown (variable) model parameters

for ct = 1:length(fixParam)
   assignin('base', fixParam{ct}, fixParamNum(ct));
 end

% Truth model response
opts = simset('InitialState', [], 'OutputPoints', 'specified', ...
	      'Solver', 'ode4', 'FixedStep', tsamp, ...
	      'SrcWorkspace', 'Current', 'Trace', '');
[T,X,Y] = sim(system, [t(1) t(end)], opts, [t,u]);

assignin('base', 'var', 0.0); % No artificaial noise in simulations
P = LocalOptimize(varParam, fixParam, fixParamNum, system, T,Y, t,u);
e = abs((exactP-P)./exactP*100);
disp(sprintf('\nParameter error: cv: %4.2f%% \t fd: %4.2f%%', e(1), e(2)));


% --------------------------------------------------------
function P = LocalOptimize(varParam, fixParam, fixParamNum, system, T,Y, t,u)
% Optimization data
X0 = [0, 0];
LB = [0, 0];
UB = [inf, inf];
options = optimset('Display', 'iter');

% Non-linear least-squares optimization
P = lsqnonlin(@LocalCostFcn, X0,LB,UB,options, varParam, T,Y,t,u, system);


% --------------------------------------------------------
function F = LocalCostFcn(varParamNum, varParam, T, Y, t, u, system)
for ct = 1:length(varParam)
  assignin('base', varParam{ct}, varParamNum(ct));
end

% Set simulation parameters.
opts = simset('InitialState', [], 'OutputPoints', 'specified', ...
	      'Solver', 'ode4', 'FixedStep', t(2)-t(1), ...
	      'SrcWorkspace', 'Current', 'Trace', '');
[T1,X1,Y1] = sim(system, [t(1) t(end)], opts, [t,u]);
F = (Y1-Y);