Equation d��tat d�un mod�le de transport-diffusion

Document Sample
Equation d��tat d�un mod�le de transport-diffusion Powered By Docstoc
					Equation d’état d’un
modèle de transport-
    diffusion.
   Applications

      Gilles Roussel,
      Eric Ternisien

    LASL (EA2600) / ULCO
           Calais
                           1
                   Plan
I   Contexte

II Sans Modèle / Avec Modèle ?

III Modélisation

IV Discrétisation, codage d’état

V Propriétés

VI Applications

VII Conclusion

                                   2
                I.          Contexte

• Problématique Inverse pour la surveillance
  de la pollution
       Mesures aux capteurs de concentration  déterminer :
                 » Le flux à la (aux) source(s) (déconvolution,
                    séparation)
                 » La position de la source (localisation)
                 » les paramètres du modèle (identification)


• Domaines applicatifs
   • Surveillance de la pollution atmosphérique
            Développer une aide au diagnostic des pics
              d’émission canalisée d'origine industrielle en
              identifiant la cause (flux, position)
   • Surveillance de la pollution fluviale ou
     phréatique
            Localiser les sources de pollution liée à des rejets
              ou fuites
   • Localisation des essais nucléaires (Traité du
     CTBTO)
            Détection des explosions  1Ktonne eqTNT à partir
              des données de concentration en xénon
            Localiser les sources de pollution liée à des rejets 3
              ou fuites
    II. Sans modèle / avec modèle ?

• Méthodes sans modèle “source - capteurs” :
    – Interpolation spatiale à partir des concentrations
      observées
         • nb capteurs ?
         • sensibilité aux paramètres ?
         • portabilité du comportement à d’autres sites ?.
    – Prédiction temporelles à partir des séries
         • déconvolution ?
         • localisation ?
         • sensibilité au paramètres ?
    Conclusion : combinaison espace - temps difficile
• Méthodes avec modèle “source - capteurs” :
    – Modèle de “comportement”
         • choix ?
         • prise en compte des paramètres physiques ?
         • portabilité ?
         • localisation?

    – Modèle de “connaissances”


                                                             4
   II. Avec modèle /sans modèle ?

• Modèle de “connaissances” :
    –   phénomènes physiques caractéristiques
    –   précision évolutive
    –   prend en compte les variables pertinentes
    –   modélisation espace-temps explicite


• Souhaits pour le modèle directe
    – simplicité, fidélité : choix ?
    – faire apparaitre explicitement la (les) source(s) et les
      variables observables
    – garder la dimension spatiale et temporelle (mesures
      spatialement réparties, cas multisources)
    – possibilité de se placer dans un cadre discret
      (mesures échantillonnées)




                                                             5
                        III. Modélisation

      • Pourquoi le modèle de transport-diffusion ?
          – Modèle (historique) de dispersion de la pollution
            dans l’air
          – bonne approximation de la propagation horizontale
          – hypothèses sur la dispersion turbulentes verticales
          – Modèle directe dynamique
          – solution en statique classique : “équation panache”


      • Le modèle de transport-diffusion (2D)
          – il satisfait en tout point la conservation de la masse :


    ( x, y, t )  ( x, y, t )      ( x, y, t )   ( x, y, t )   ( x, y, t )
                 Ux             Uy                K               K
        t            x                 y         x    x         y    y
                              advection                  diffusion

  (x,y,t)        : concentration au point (x,y) à l ’instant t
U  [Ux,Uy]T      : vent

  K               : diffusion


                                                                                        6
               III. Modélisation

• Hypothèses
   – Domaine atmosphérique non borné
   – Réduction à un domaine de calcul (de surveillance)
      borné de contour artificiel 
   – conditions aux limites

                                                          Atmosphère
                              Ex. conditions de Neuman :
                              continuité aux limites artificielles
                     n
                                                 
                


                     Source
                                     ’

              Vent                Ex. conditions de Dirichlet



                     Vue de dessus


              Vent
                              
      Terre                                                          7
                              III. Modélisation

• Conditions aux limites (locales)
      – Limites artificielles non réfléchissantes
         • continuité (Pearson)
                          (x,y,t) (x,y,t) (x,y,t)                                                          Sur 
                                   x         y        0
                             t         x          y
                           – x,  y célérité équivalente
      – Limites physiques réfléchissantes
                 • exemples =0 sur ’


      – à une source : au point (xs,ys)

    ( xs , ys , t )      ( xs , ys , t )       ( xs , ys , t )   ( xs , ys , t )   ( xs , ys , t )
                      Ux                    U y                    K                   K                    S ( xs , y s , t )
         t                     x                      y           x      x           y      y



• Conditions initiales : pollution de fond
                                                                                                       (0, x, y )  f ( x, y )

• Modèle elliptique linéaire




                                                                                                                                 8
                 IV. Discrétisation
• On discrétise les signaux S(xs,ys,t) et  (x,y,t) avec
  le même pas d'échantillonnage temporel T
• On discrétise l'espace avec les pas d'échantillonnage
  directionnels h et k
        Titre:
        fig11.eps
        Auteur:
        fig2dev Vers ion 3.2 Patchlevel 1
        Aperçu:
        Cette im age EPS n'a pas été enregistrée
        avec un aperçu intégré.
        Commentaires :
        Cette im age EPS peut être im primée sur une
        imprimante PostScript mais pas sur
        un autre type d'im primante.




• Les méthodes de discrétisation:
    – différences finies,
    – élements finis


                                                       9
               IV. Discrétisation

• Dérivation par les différences finies en un point i:


                                                     
      O                   i-1n                                    X
                                           in             i+1n


    – dérivée première temporelle :
                                                 n1
            i  (t  T )   (t )           i     n
                                    0(T )  i
            t          T                       T
    – dérivée première spatiale
          • centrée espace, schéma de Lax : effet de lissage
          • décentré en espace (“upwind”):
                              in i1n
                            h siU.x0
                            
                          x  i1n in
                                         siU.x0
                              h

                – calcul effectué que pour les points "au vent”
                – tenir compte du sens du vent
                – Cette méthode ajoute un terme de viscosité négligeable si
                  le nombre de Reynolds de maille Rh:


                                     U .h
                              Rh           1 5
                                      K


                                                                       10
              IV. Discrétisation

– Dérivée spatiale seconde
     • explicite : in+1=f1( in,  vn),  v potentiel au voisinage de i
             – Markovien,

   2 i  (t , x  h)  2 (t , x)   (t , x  h)             2 i   i 1
                                                               n       n      n
                                                    0(h2 )  i 1
   x
    2
                            h 2
                                                                    h2


     • implicite :  in+1 =f2( in ,  in+1,  vn+1)
             – schéma stabilisant,
             – résolution numérique supplémentaire à chaque n
             – non Markovien


     • Crank-Nicolson :  in+1 =f1(.) + (1-  )f2(.)
             – idem implicite




                                                                                  11
                                           IV. Discrétisation

        • Equation d'advection - diffusion (cas 2D,
          K=cste)
       in1  m1in 1, j  m2in j  m3in 1, j  m4in j 1  m5in j 1
         ,j                    ,                     ,            ,                                       (i, j)  intérieur ( )


        • conditions limites
                      – aux frontières artificielles (condition de continuité)

       in1  m7in j  m8in 1, j  m9in j 1
         ,j        ,                     ,                             (i, j)  

                      – aux sources
        in1  m1in 1, j  m2inj  m3in 1, j  m4inj 1  m5inj 1  m6 Sin, j
          ,j                    ,                    ,           ,                                             en (i  i s , j  js )


                      – obstacle (sans turbulence)

                             in1  0
                               ,j               (i, j)   ' '


        • Paramètres (pour Ux, Uy positifs)                                                                y
                                                                                           m7  1  T  x   
                                                                                                       h     k 
    TK
m1  2 x                                                                TK y                                   
     h
                      TU y
                                                                 m4                           1  (m8  m9 )  m2  
           TU x                  2TK x 2TK y          5
                                                                         k2
m2  1                                    1   mi                                             x
            h          k          h2    k2        i 1,i  2            T Uy       TK y    m8  T             m3  
       TU x TK x                                                 m5                                   h
m3         2                                                            k         k2             y
        h    h                                                   m6  T                   m9  T             m5  
                                                                                                    k

        • Précision
                                  •       e=O(T)+O(h2)+O(k2)+ O(h)+O(k)
                                                                                                                           12
                      IV. Discrétisation

                                                      capteurs
                J=L
Contour artificiel                  c1
                                                                 Obstacle ’
non réfléchissant                                c2
                     c3                i,j+1
                                        *
                            i-1,j    * i,j * i+1,j
                                        *
                                      i,j-1
      source

               i=1                                    i=l

 • Schéma numérique
        – en 




        – en              0      0     0
                            m      m7    0
                             8            
                            0
                                   m9    0
                                           




                                                                               13
                                             IV. Discrétisation
           • Modèle d’état
                                  vecteur d’état                                                                                                         n T
                                                                            n  [1,1 ,......, 1, ],[2,1 ,......, 2, ],........ .,[L,1 ,......, L, ]
                                                                                        n            n        n            n                 n
                                                                                                                                                            
                                  Equation d’état                           n1  A n  T .BS n
                                                                            n
                                                                            Y  C n
                                                                                                                                      m4               
   A1,1       A1, 2                                                            m2        m1                         Ak ,k 1     
                                                                                                                                                      
                                                                                                                                                        
  A                                                                            m                                                                m4 
   2,1        A2, 2   A2,3                                                      3        m2     m1                                                 
                                                                   
A
                               
                                                                       Ak ,k                                    
                                                                                                                            m5                  
                                                                                                                                                      
                                                          AP 1, P                              m3 m2           m1 
                              AP 1, P  2   AP 1, P 1                                                                A                          
  
                                             AP , P 1    AP , P 
                                                                                
                                                                                                     m3               k 1,k 
                                                                                                                   m2                                  
                                                                                                                                
                                                                                                                                                    m5 
                                                                                                                                                        
  A caractérise l’évolution libre du système de dispersion
        0 
                                                                                                             0  0 1 0   0 
                                                                                                            C  0 1 0     0 
               B modélise l’implantation des sources
        0                                                                                                                     
                                                                                                              0    0 1 0 0 
    B  1                                                                                                                     
        0 
                 dim(B)= lL*ns
         
        
                                                                                C modélise l ’implantation des capteurs
        0 
         
           • Les bruits                                                         dim(C)= nc*lL
                       – Wn : bruit d’état
                                  • modélise l’incertitude sur l’état (turbulence non prise
                                    en compte par le modèle, topographie)
                                  • moyenne parfois non nulle, trace(E(Wn.WnT)) faible
                                  • indépendant de Sn et Vn

                       – Vn : bruit de mesure
                                                                                                           n1  A n  T .BS n  wn
                                  • moyenne souvent nulle                                                  n
                                                                                                           Y  C  v
                                                                                                                        n     n
                                  • trace(E(Vn.VnT)) faible
                                                                                                                                                     14
                    IV. Discrétisation

• Stabilité du schéma numérique
   – condition nécessaire de convergence de la solution
     i,jn.
   – peut s’étudier à partir des valeurs propres de la
     matrice A
   – rayon spectral (A)=max|k| 1kLl+1
   – on peut montrer (Gershgörin-Hadamard)
              L.l
      k     ai , j  1  k  lL
              j 1

                       L.l           5
   – comme             ai , j     mk  1      (A matrice de Markov)
                       j 1       k 1



   – alors        ( A)  max  k  1                 1  k  lL

   – finalement, on peut montrer

                         ( A)  1       si m2  0



   d'où la condition sur la période d’échantillonnage T
      pour un pas spatial (h,k) donné
                                           h2k 2
            ( A)  1  T  2
                           2k K x  2h2 K y  hk 2 U x  kh 2 U y



                                                                         15
• Simulations
      Titre:
                      ig17.eps
      T:\latex\Thes e\f
      Auteur:
      MATLAB, The Mathworks , Inc.
      Aperçu:
      Cette im age EPS n'a pas été enregistrée
      avec un aperçu intégré.
      Commentaires:
      Cette im age EPS peut être im primée sur une
      imprimante PostS cript mais pas sur
      un autre type d'im primante.




                                                     sans advection
  Titre:
  T:\latex\Thes e\fig16.eps
  Auteur:
  MATLAB, The Mathworks , Inc.
  Aperçu:
  Cette im age EPS n'a pas été enregistrée
  avec un aperçu intégré.
  Commentaires:
  Cette im age EPS peut être im primée sur une
  imprimante PostScript mais pas sur
  un autre type d'im primante.




                                                     avec advection
                                                                      16
                         V. Propriétés

• Observabilité
      – Observabilité stricte, c’est la possibilité de résoudre :
                    y ( k0 )   C 
                    y (k  T )   CA 
                         0            (k )  Y  .(k0 )
                         ...      ...  0
                                  n1 
                    y (k0  nT ) CA 

      – O doit être de rang = l.L = dimension de 
      – un nombre de mesures égale à :
            • cas monocapteur : dim(Y)  l.L
            • cas nc capteurs : dim(Y)  (l.L/nc)


      – si (rang(O)=r )< l.L seulement (l.L-r) noeuds ne
        peuvent être estimés.
                                                                   Titre:



      – Observabilité théorique:
                                                                   fig15.eps
                                                                   Auteur:
                                                                   fig2dev Version 3.2 Patchlevel 0-beta3
                                                                   Aperçu:
                                                                   Cette im age EPS n'a pas été enregis trée
                                                                   avec un aperçu intégré.
                                                                   Commentaires :
                                                                   Cette im age EPS peut être im primée s ur une

                  nclL
                       , lL)  demi largeur de bande de A
                                                                   imprimante PostScript mais pas sur


rang (O)  inf(
                                                                   un autre type d'im primante.



                   2


      – Le rang est meilleur dans le cas multicapteurs
      – numériquement : rang(O)  card (i  eig (O) / i  seuil)

      – le rang renseigne sur la qualité de restauration de 
              (k0 )  O Y ,    O  inverse au sens de Moore Penrose
                                                                                                17
                 Titre:
                 H:\LATEX\These\Fig192.eps
                 Auteur:
                 MATLAB, The Mathworks , Inc.
                 Aperçu:
                 Cette im age EPS n'a pas été enregistrée
                 avec un aperçu intégré.
                 Commentaires :
                 Cette im age EPS peut être im primée sur une
                 imprimante PostScript mais pas sur
                 un autre type d'im primante.




•     Evolution de rang(O) pour 1 capteur, seuil de rang :1e-25

            Titre:
            H:\LATEX\These\Fig193.eps
            Auteur:
            MATLAB, The Mathworks , Inc.
            Aperçu:
            Cette im age EPS n'a pas été enregis trée
            avec un aperçu intégré.
            Commentaires :
            Cette im age EPS peut être im primée s ur une
            imprimante PostScript mais pas sur
            un autre type d'im primante.




    Evolution de rang(O) pour 2 capteurs, dont1 fixe en (8,8), seuil de
                                rang :1e-25
•     Conclusion :
        – plus le capteur est éloigné de la source dans la direction du
          vent, meilleur est l'observabilité
                                                                      18
        – plus il y a de capteurs, meilleur est l'observabilité
                                                         V. Propriétés
• Commandabilité (Excitabilité des nœuds)

   – Existe t-il une commande S  s(k0 ), s(k1 ),........ .....s(k f  1)T

                                                         vérifiant :
                                                 s ( k0 ) 
                                                 s(k ) 
             X (k f )  A B,......... ..., AB B 
                         n 1
                                                
                                                       1     C.S
                                                            
                                                                       
                                                           
                                                 s(k f  1)

   La réponse est oui si rang(C) = l.L, sinon il y a (rang(C)-l.L)
      nœuds non excités.
   – On peut montrer:
                                                           nslL
     rang (C )  inf(                                           , lL)  demi largeur de bande de A
                                                            2
   – De façon logique, le rang(C) augmente avec le nombre de
     source.
   – Un état X(kf) peut être atteint en n coups (n<l.L), si
         •            X(kf) appartient au sous espace engendré par les colonnes de

                                                                             
                      C
                                 C  An1B,......... ..., AB B
         • La représentation des colonnes de C sur le maillage, constitue
           une base de motifs
          Titre:
          fig195.eps
          Auteur:
          fig2dev Vers ion 3.2 Patchlevel 1
          Aperçu:
          Cette im age EPS n'a pas été enregistrée
          avec un aperçu intégré.
          Commentaires:
          Cette im age EPS peut être im primée sur une
          imprimante PostScript mais pas sur
          un autre type d'im primante.




                                                                                                19
                                      V. Propriétés
• Commandabilité (suite)
       –   En pratique, le rang(C) dépend du seuil de sélection des valeurs
           propres considérées non nulles.
       –   Pour un seuil d'excitabilité donné, le rang indique le nombre de
           nœuds excitables.


              Titre:
              H:\LATEX\These\Fig197.eps
              Auteur:
              MATLAB, The Mathworks , Inc.
              Aperçu:
              Cette im age EPS n'a pas été enregistrée
              avec un aperçu intégré.
              Commentaires :
              Cette im age EPS peut être im primée sur une
              imprimante PostScript mais pas sur
              un autre type d'im primante.




    Evolution de rang(C) en fonction de la position de S, à seuil fixé

•    Conclusion : Pour une direction de vent donnée, plus la
     source est au vent, plus il y a d'état excités.




                                                                              20
                    V. Propriétés

• Identifiabilité
    – C.N. : Commandabilité et Observabilité



    on augmente la probabilité d'identifier le modèle si :
         • le capteur est placé sous le vent, le plus loin de la source
         • le nombre de capteurs augmente

                                Ouf !!!




                                                                          21
                   VI Applications

   • Principe Majeur :
         – Considérer la relation source-capteur comme un filtre de
           convolution

                                      v1(n)     Bruit
                                                capteur
                            H1 (z)
                                                Y1(n)
                                         


  Source                                                      mesures
  S(n)
   Position

                           H2 (z)              Y2(n)
                                         

                                     v2(n)




Hi(z)=C(zI-A)-1B = f [position(S), position(capteursi), K, U]            i




                                                                        22
                      VI Applications
• Différentes formes du modèle :
                                                                 = paramètres physiques
    – forme d'état A(1), B ( 2), C 1
                                     2 = paramètre de localisation
    – réponses impulsionnelles

   h( j ) (1, 2 )  [h0 (1, 2 ), h1 (1, 2 ),........ .....hM          (1, 2 )]
                       ( j)            ( j)                          ( j)



                    hi (1, 2 )  CAi 1 (1 ) B( 2 )
                      ( j)
                                                                                         M = ordre RIF
    – matrice de Hankel
                      H (1,2 )  O(1).C(1,2 )
    – matrice de filtrage
                                                    h0        hM 0  0 
                                                                         
                                                   0                
                                              HN  
                                                                   0 
                                                                         
                                                   0          0 h0  hM 
                                                                         

• Différents problèmes inverses :
    – Estimation des paramètres physiques 1
    – Localisation de source(s) : estimation de 2 (i.e. B)
    – Déconvolution / séparation de sources : estimation de
      Sn
    – Prédiction temporelle
• Informations :
    – Mesures capteurs, estimation de bruit, entrées
    – Mesures capteurs, estimation de bruit, info a priori
                                                                                                  23
                            VI Applications




Exemple :
        Estimation des réponses impulsionnelles
        par méthode spectrale paramétrique
        (séparation des sous-espaces source et bruit)



 Titre:
 I:\latex\These\fig29-4.eps
 Auteur:
 MATLAB, The Mathworks , Inc.
 Aperçu:
 Cette im age EPS n'a pas été enregistrée
 avec un aperçu intégré.
 Commentaires :
 Cette im age EPS peut être im primée sur une
 imprimante PostScript mais pas sur
 un autre type d'im primante.




                                                        24
             VII     Conclusion


• Modélisation d'un phénomène de transport-
  diffusion par une représentation d'état.

• réutilisabilité de la méthodologie ?

• Point de départ pour bon nombre de problèmes
  d'estimation

• Aspects applicatifs : surveillance de points sources




                                                     25

				
DOCUMENT INFO
Shared By:
Categories:
Tags:
Stats:
views:3
posted:11/26/2011
language:French
pages:25