Thread Subject:

Using pem with a non smooth dynamical system

Hi,

I'm stuck with parameter estimation of an idnlgrey model. A similar model is the bouncing ball, with gravity g and damping coefficient c. How can I define an odefile with nonterminal events ? User-defined resets failed with stiff solvers. I tried to reuse outputs from the Stateflow demo sf_bounce, but parameters were not updated:
Warning: Covariance matrix estimate unreliable. Not stored.

Any ideas ?
Cheers

Julien

 ----
global ys dxs ts

load_system('sf_bounce');
set_param('sf_bounce/C1', 'Value', strcat('[10 0.1]'));
set_param('sf_bounce/C2', 'Value',num2str(0.8));
set_param('sf_bounce/C3', 'Value',num2str(-9.81));
set_param('sf_bounce/C4', 'Value',num2str(9.81));
sim('sf_bounce');
t0=ysim.time;
y1=ysim.signals.values;

ys=[];
dxs=[];
ts=[];

z = iddata(ysim.signals.values, [], 0.01, 'Name', 'SCV');
set(z, 'OutputName', 'Position');
set(z, 'OutputUnit', 'cm');
set(z, 'Tstart', 0, 'TimeUnit', 's');

% non-smooth model

FileName = 'ball_m';
Order = [1 0 2];% Model orders [ny nu nx].
Parameters = [8.5 0.75];%g=9.81 c=0.8
InitialStates = [10; 0.1];% Initial parameter vector.
Ts = 0.01; % Initial initial states.
nlgr = idnlgrey(FileName, Order, Parameters, InitialStates, 0);

set(nlgr, 'OutputName', 'Position', ...
          'OutputUnit', 'cm' ,...
          'TimeUnit', 's');
setinit(nlgr, 'Name', {'Position' 'Speed'});
setpar(nlgr, 'Minimum', {5 0});
setpar(nlgr, 'Maximum', {10 1});

figure
compare(z,nlgr,3);
title('After')

nlgr = pem(z, nlgr, 'Display', 'Full','MaxIter',100); % Perform parameter estimation.

figure
compare(z,nlgr,3);
title('After')

disp(nlgr.Parameters(1).Value)
disp(nlgr.Parameters(2).Value)

function [dx, y] = ball_m(t,x,u,g,c,varargin)
    global ys dxs ts

        if isempty(ys)
            disp('test')
            y0=[10,0.1];

            load_system('sf_bounce');
            set_param('sf_bounce/C1', 'Value', strcat('[',num2str(y0),']'));
            set_param('sf_bounce/C2', 'Value',num2str(c));
            set_param('sf_bounce/C3', 'Value',num2str(-g));
            set_param('sf_bounce/C4', 'Value',num2str(g));
            sim('sf_bounce');
            ts = ysim.time;
            ys = ysim.signals.values;
            dxs=dxsim.signals.values;
        end
        
        y = interp1(ts,ys,t);
        dx = interp1(ts,dxs,t)';

 end