Transcript Commission de discipline - Ligue d`Alsace de Handball

TP CONTSID
Hugues Garnier
Estimation param´
etrique de mod`
eles `
a temps continu
`
a partir de donn´
ees ´
echantillonn´
ees
Initiation `
a la boˆıte `
a outils CONTSID pour Matlab
A rendre le vendredi 24 janvier 2014
1
Etude en simulation
Identification d’un mod`ele `
a temps continu du 2i`eme ordre par LSSVF, IVSVF et SRIVC
1.1
G´
en´
eration de signaux d’entr´
ee/sortie
Soit le syst`eme d´ecrit par :
S : y(tk ) =
B(p)
u(tk ) + e(tk )
F (p)
(1)
o`
u p repr´esente l’op´erateur de diff´erentiation, u l’entr´ee du syst`eme, y la sortie mesur´ee, e un bruit
blanc `a temps discret, et les polynˆ
omes B(p) et F (p) sont d´efinis par :
(
B(p) = 2
F (p) = p2 + 4p + 3
1. D´eduire de l’´equation (1) la structure du mod`ele du syst`eme (CARX ou COE)
2. Donner les caract´eristiques du syst`eme `a temps continu : gain statique, constantes de temps
et pulsation de coupure associ´ee `
a chaque constante de temps
3. Pr´eciser l’´equation diff´erentielle du syst`eme `a l’instant t = tk
Vous allez dans un premier temps g´en´erer des signaux d’entr´ee/sortie `a partir desquels vous allez
tenter d’estimer les param`etres du syst`eme (1). Pour cela, une fois sous Matlab, entrer, dans un
script Matlab, la s´equence de commandes suivante :
N=1000;
B=[2];
F=[1 4 3];
Te=0.01;
t=(0:N-1)’*Te;
u=square(2*pi*0.5*t);
randn(’state’,sum(100*clock))
e=0.2*randn(N,1);
x=lsim(B,F,u,t);
y=x+e;
datadet=iddata(x,u,Te);
data=iddata(y,u,Te);
1
figure(1)
idplot(datadet)
figure(2)
idplot(data)
Vous avez `a pr´esent `
a votre disposition deux jeux de donn´ees contenus dans les objets de donn´ees
datadet et data pour effectuer l’estimation param´etrique.
1.2
1.2.1
Estimation param´
etrique
Choix de la structure de mod`
ele
On supposera les ordres du mod`ele connus.
1.2.2
Estimateur LSSVF
On s’int´eresse dans un premier temps `
a l’identification du syst`eme par la m´ethode des moindres
carr´es/filtre de variable d’´etat (LSSVF). Le mod`ele recherch´e est du type ARX continu :
MCARX : A(p)y(tk ) = B(p)u(tk ) + e(tk )
(2)
1. Rappeler l’´equation diff´erentielle associ´ee au mod`ele.
2. En d´eduire l’´ecriture du mod`ele sous forme de r´egression lin´eaire `a l’instant t = tk .
3. En supposant N mesures des signaux d’entr´ee/sortie, formuler le probl`eme sous forme matricielle. Donner la dimension de chaque matrice.
4. Rappeler la solution de l’estimateur LSSVF (voir transparents de cours).
5. Implanter cette solution sous Matlab dans le cas des donn´ees d´eterministes puis bruit´ees. Commentaires. On rappelle que l’on r´ealise ici une estimation ponctuelle et que l’on ne peut parler
de biais de l’estimateur mais plutˆ
ot d’erreur d’estimation. Vous pouvez ex´ecuter plusieurs fois
votre programme pour d´egager la tendance au niveau du biais des estim´ees.
6. La boˆıte `a outils CONTSID 1 contient diff´erentes fonctions permettant de r´ealiser l’estimation
des param`etres d’un mod`ele. Elles ont toutes la structure suivante :
M=fonction(data,nbreparam,hyperparametre)
avec :
• M : le mod`ele identifi´e ;
• fonction : le nom de la fonction CONTSID : lssvf, ivsvf, srivc, . . .
• data : l’objet donn´ees contenant les signaux d’entr´ee/sortie construit `a la prtir de la
fonction iddata : data=iddata(y,u,Te) ;
• nbreparam : un vecteur ligne permettant de d´efinir le nombre de param`etres des
polynˆ
omes du mod`ele et la valeur du retard (suppos´e ˆetre un multiple entier de la
p´eriode d’´echantillonnage τ =nk×Te ) : [na nb nk] par exemple. Attention il s’agit du
nombre de param`etres `
a estimer de chaque polynˆome du mod`ele et non du degr´e de
chaque polynˆ
ome.
• hyperparametre : hyper-param`etre ´eventuel de la m´ethode d’estimation : pulsation de
coupure (rad/s) du filtre SVF par exemple.
1
www.cran.uhp-nancy.fr/contsid/
2
La commande suivante permet d’afficher les param`etres du mod`ele estim´e :
present(M);
Si vous d´esirez obtenir plus d’informations sur l’utilisation d’une fonction, vous pouvez entrer
la commande suivante `
a tout moment :
help fonction
La fonction CONTSID lssvf permet d’estimer les param`etres d’un mod`ele `a temps continu
par LSSVF. La commande est la suivante :
Mlssvf=lssvf(data,[na nb nk],lambda);
na et nb correspondent respectivement au nombre de param`etres `a estimer des polynˆomes
A(p) et B(p) et nk repr´esente le nombre d’´echantillons pour le retard, lambda repr´esente la
pulsation de coupure du filtre SVF (en rad/s).
Les param`etres estim´es peuvent ˆetre visualis´es par la commande present(Mlsssvf) :
7. Entrer les commandes suivantes pour estimer les param`etres du mod`ele par la m´ethode des
moindres carr´es/filtre de variable d’´etat (LSSVF) :
Mlssvf1=lssvf(datadet,[2 1 0],3);present(Mlssvf1);
Mlssvf2= lssvf(data,[2 1 0],3);present(Mlssvf2);
8. Pour une r´ealisation du bruit, comparer les r´esultats de votre implantation avec ceux obtenus
via la fonction de la boˆıte `
a outils CONTSID.
1.2.3
Estimateur IVSVF
On s’int´eresse `
a pr´esent `
a l’identification d’un mod`ele CARX par la m´ethode de la variable
instrumentale `a mod`ele auxiliaire/filtre de variable d’´etat IVSVF.
1. En supposant N mesures des signaux d’entr´ee/sortie. Rappeler la solution de l’estimateur
IVSVF. Pr´eciser la dimension de chaque matrice.
2. Implanter cette solution sous Matlab (voir transparents de cours). Commentaires.
3. La fonction ivsvf permet d’estimer les param`etres d’un mod`ele CARX par la m´ethode de la
variable instrumentale associ´ee `
a la technique des SVF. La commande est la suivante :
Mivsvf=ivsvf(data,[na nb nk],lambda);
4. Entrer les commandes suivantes pour estimer les param`etres du mod`ele par IVSVF :
Mivsvf= ivsvf(data,[2 1 0],3);
present(Mivsvf);
5. Comparer avec les r´esultats obtenus avec lssvf sur un mˆeme jeu de donn´ees.
3
1.3
Estimateur SRIVC
La fonction srivc permet d’estimer les param`etres d’un mod`ele du type COE :
MCOE : y(tk ) =
B(p)
u(tk ) + e(tk )
F (p)
(3)
par la m´ethode it´erative de variable instrumentale. Les estim´ees d´elivr´ees sont alors optimales (sans
biais et `a minimum de variance) lorsque le bruit additif e(tk ) sur la sortie est blanc (syst`eme vrai
de type COE). La commande est la suivante :
Msrivc=srivc(data,[nb nf nk]);
% Attention `
a la d´
efinition du nombre de param`
etres `
a estimer
nb et nf correspondent respectivement au nombre de param`etres `a estimer des polynˆomes B(p) et
F (p) du mod`ele et nk repr´esente le nombre d’´echantillons pour le retard. Les param`etres estim´es
ainsi que leur ´ecart-type peuvent ˆetre visualis´es par la commande :
present(Msrivc);
1. Entrer les commandes suivantes pour estimer les param`etres du mod`ele COE `a l’aide de la
fonction srivc :
Msrivc=srivc(data,[1
2
0]);present(Msrivc);
2. Comparer avec les r´esultats obtenus avec lssvf et ivsvf.
1.3.1
Analyse comparative des performances des estimateurs LSSVF, IVSVF et
SRIVC par simulation de Monte Carlo
Ecrire un programme qui permet de comparer les performances des estimateurs LSSVF, IVSVF et
SRIVC par simulation de Monte Carlo dans le cas du syst`eme (1).
1.3.2
Programmes de d´
emonstration de la boˆıte `
a outils CONTSID
Pour acc´eder aux programmes de d´emonstration disponibles dans la boˆıte `a outils CONTSID, entrer
idcdemo
Ex´ecuter les diff´erents programmes de d´emonstration disponibles.
2
Applications sur des donn´
ees r´
eelles
Cette partie a pour objectif de mettre en application la m´ethodologie de l’identification de mod`eles
`a temps continu sur des donn´ees r´eelles.
2.1
Identification d’un syst`
eme thermique
1. Protocole exp´erimental - Les donn´ees brutes ´echantillonn´ees `a Te=80ms se trouvent dans le
fichier dryer2.mat. Celles-ci peuvent ˆetre charg´ees en entrant la commande suivante :
load dryer2;
On choisit d’utiliser les 500 premiers points pour estimer un mod`ele. On forme l’objet donn´ees
(data) contenant les donn´ees d’entr´ee/sortie brutes :
data=iddata(y2(1:500),u2(1:500),0.08);
4
Tracer l’´evolution temporelle des signaux d’entr´ee/sortie :
idplot(data);
2. Pr´e-traitement des donn´ees - Dans un premier temps, il faut ´eliminer les composantes continues des signaux :
data=detrend(data,’const’);
idplot(data);
3. Choix de la structure et d´etermination des ordres/retard du mod`ele - En utilisant la proc´edure
de recherche de structure de mod`eles (voir idcdemo8), proposer une structure et s´electionner
les ordres/retard du mod`ele
4. Estimation param´etrique - D´eterminer les param`etres du mod`ele choisi par la m´ethode de
votre choix.
5. Validation - Nous allons juger de la qualit´e du mod`ele obtenu en consid´erant une portion
de la campagne de donn´ees qui n’a pas ´et´e utilis´ee pour r´ealiser l’identification (validation
crois´ee).
uv=detrend(u2(501:900),’const’);
yv=detrend(y2(501:900),’const’);
dataval=iddata(yv,uv,0.08);
Afin de v´erifier si le mod`ele d´ecrit bien le syst`eme, on effectue une simulation du mod`ele puis
on compare les sorties du syst`eme et du mod`ele :
comparec(dataval,M);
V´erifier l’ad´equation du mod`ele.
6. Poursuivez la validation en effectuant les autres tests habituels (analyse de corr´elation des
r´esidus, trac´e des pˆ
oles et des z´eros).
2.2
Identification d’un m´
ecanisme pi´
ezo´
electrique
Les donn´ees ´echantillonn´ees `
a Te=33 × 10−6 s se trouvent dans le fichier piezo.mat.
1. Appliquer la proc´edure d’identification pour d´eterminer un mod`ele mod`ele `a temps continu
en utilisant une approche directe ou indirecte du m´ecanisme pi´ezo´electrique.
2. Justifier et expliquer tous vos choix.
3. Donner les caract´eristiques principales du mod`ele que vous proposez : ordre, gain statique,
pulsation(s) propre(s) non amortie(s), coefficient(s) d’amortissement(s).
5