Discussion :

[Matersss]

Bonjour,

J'utilise des matrices très mal conditionnées A1 et A2 (qui sont symétriques définies positives).
Je dois prendre les valeurs propres ce ces produits de matrices :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
[V1,D1]=eig(A2\A1)
[V2,D2]=eig(A1\A2);

Le problème est que les résultats sont en général très mauvais et j'ai les warnings de Matlab pour me dire que mes matrices sont mal conditionnées (avec des rcond atteignant 10^-20).
Si au lieu d'utiliser

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

eig(A1\A2)

, j'utilise

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

[V2,D2]=eig(inv(A1)*A2);

, c'est un peu mieux.
Sinon, j'ai déjà essayé d'utiliser

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

[V2,D2]=eig(A2,A1);

vu que ça revient à un problème de valeurs propres généralisé.
Je cherche donc des moyens d'améliorer le conditionnement de mes matrices sur Matlab (comme "balance" mais cela ne marche pas pour les matrices symétriques définies positives) ou alors de calculer plus précisément les inverses et valeurs propres de mes matrices.

Merci
Matersss

[FR119492]

Salut!
Est-ce que le mauvais conditionnement de tes matrices provient de la nature du problème que tu dois résoudre ou plutôt d'une mauvaise formulation de ce problème?
Jean-Marc Blanc

[tachmou]

Bonjour,

Je pense que tu touches du doigt à des problèmes vachement complexes... sache que la recherche de valeurs propres de matrices mal conditionnées... est encore un problème pas totalement résolu... et pour le peu que j'ai entendu dire : pour ce qui est résolu, ça peut être très gourmand en temps de calcul, surtout si tu n'utilises pas de noyau assembleur (et matlab ne l'a pas implémenté, même pas dans les dernières versions)... En tout cas, perso, je ne peux pas t'aider à fond...

Cependant, petites questions de curiosité :
- tes matrices sont bien carrées ? (si non : utiliser plutôt PINV et renseigne toi sur ce que ça fait)
- taille de tes matrices ?
- tu travailles sur quel problème ? (quelle est le degré de physique impliqué)
- aurais-tu à tout hasard une idée de si tu peux transformer tes matrices sous forme diagonale par blocs ? (ou tout du moins avoir les éléments de gros poids sur la diagonale et/ou autour de la diagonale : si j'ai bien compris, ça peut aider )

Cordialement,

[shaiHulud]

Bonjour,

Quelques pistes:

* ACP
Forcer l'isotropie sur les petits espace propres, en remplaçant B'*D*B par B'*Di*B; ou je force les K plus petites valeurs propres de Di à être égales, tout en conservant Trace(D)= Trace(Di). L'inversion devient trivial sur l'espace propre isotrope.

* Réduction d'un effet taille
Si S=var(X) (ta matrice s’interprète comme une matrice de variance d'un vecteur X) et si le 1er vecteur propre est dominant et à toutes ses coordonnée de même signe (on appelle CA un effet taille en ACP)
- Considérer les variables en différences 2 à 2 : les (X_i - Y_j)j>i plutôt que les Xi. Cette transformation amoindrit la force de la 1ere valeur propre puisque l'on se place dans l'orthogonal de (1,...,1)
- Si les variables sont positives, passer par les (Xi/Xj) j>i
- Si les variables sont positives, factoriser via une NMF (doc nnmf)

Dans les 2 derniers cas, cela t'oblige à revisiter la théorie classique basée sur des matrices de covariances, en remplaçant ta fonction de perte par une autre quantité sur laquelle les problèmes numériques sont gérables.

* Autres idées utilisées en régression linéaire (le but de l'inverse est d'inverser une relation linéaire via une projection orthogonale)
- ridge regression : inverser (S+k*Id) avec k petit. Cela biaise le résultat mais réduit sa variabilité.
- regression quantiles : inverser une sous matrice de plein rang X(h,: ) plutôt que X'*X

[FR119492]

Salut!

tu travailles sur quel problème ? (quelle est le degré de physique impliqué)

La question posée par tachmou est essentielle. Dans ces forums, toutes les discussions devraient être initiées par une réponse complète à cette question.
Jean-Marc Blanc

[Matersss]

Bonjour,

Merci pour toutes ces réponses.
Je cherche à concevoir des lois de commande pour un système multivariable (2 entrées / 2 sorties).
Une approche est de considérer que l'on commande deux systèmes monovariables mais cette approche fonctionne si les termes non diagonaux sont faibles par rapport aux termes diagonaux.
Dans un premier temps, le but est donc de rendre la matrice le plus possible diagonale dominante par colonne en ajoutant une matrice de pré-compensation. Une des méthodes pour trouver cette matrice s'appelle la pseudo diagonalisation (dans mon cas dynamique).
Je ne vais pas m'attarder sur la théorie mais le but est d'optimiser les rapports des termes diagonaux sur les termes non diagonaux (par colonne), et cela peut revenir à un problème de valeurs propres généralisées (afin d'éviter l'optimisation qui peut prendre du temps ou tomber sur des minimums locaux).
Le problème est que les matrices comportent des termes qui n'ont pas forcément le même ordre de grandeur (et parfois les ordres de grandeurs sont similaires sans pour autant que les matrices soient bien conditionnées).

Mes matrices A1 et A2 sont symétriques définies positives (donc carrées).
Elles peuvent être de taille 2*2, 4*4, 6*6... en fonction de l'ordre de la matrice de précompensation que j'impose.
Je commence à rencontrer des problèmes dès le cas 4*4.

Je vous remercie pour l'intérêt que vous portez à ce problème.
Je vais commencer à regarder vos réponses plus en détail !

Matersss

[Matersss]

Pour répondre à VV33D, je n'ai pas l'impression que la première méthode soit efficace (après, je ne connaissais pas et je n'ai pas trop regardé la théorie).
Mon problème ne respecte pas les hypothèses pour la deuxième méthode.
En revanche, il semble que la régularisation Tychonoff (ridge regression) que tu as proposé peut m'aider dans certains cas. Les résultats que j'obtiens sont parfois meilleurs et parfois moins bons mais ça peut faire avancer les choses dans certains cas.
(sinon, je n'ai pas trop compris ce qu'il fallait faire pour la régression quantile)

Matersss

[shaiHulud]

*ACP
Il est facile d'inverser une matrice une fois diagonalisée (i.e B' *inv(D)*B). Supposons D diagonale ordonnée décroissant. Si lmax / lmin >> 1 (rapport des valeurs propres trop fort), les incertitudes (statistique ou numériques) sur les derniers vecteurs propres seront amplifiées via le produit par inv(D).

L'idée derrière l'ACP est la suivante. Toutes les valeurs propres non significativement différentes des autres seront considérées comme égales. On se retrouve alors avec un espace propre E de dim>1 sur lequel la matrice est scalaire, et donc la base de diago dans cet espace propre est arbitraire ! On amoindrit donc les incertitudes provenant des derniers vecteurs propres dans B, car sur E, la résultat ne dépendra plus de la base de diago, qui est arbitraire.

Dit autrement, on a contraint les degrés de liberté de inv(X) pour le stabiliser. Une autre approximation serait de la forme X= B'*D*B + u*Identité, ou D n'est PAS de plein rang.

* Pour les secondes méthodes, j'ai écrit qu'il fallait que ta matrice s’interprète comme une matrice de variance d'un vecteur X, mais il suffirait qu'elle s’interprète comme une matrice de produit scalaires (S(i,j)= <Xi,Xj>) pour que ca marche.

* Les regressions quantiles sont à comparer aux moindres carrés en régression. Pour estimer le modèle Y=X*b+u (u indépendant de X et iid), les MCO proviennent du programme min_b E[(Y-X*b)^2] qui donne b*= inv(X'*X)*X'*y. Les régressions quantiles estiment la médiane conditionnelle et non l'espérance conditionnelle. Le programme est min_b E |Y-X*b|, et la solution est de la forme b*= inv(X(h,: ))* X(h,: )*y(h), ou h sont des indices permettant d'extraire une matrice carrée inversible X(h,: ) depuis la matrice rectangulaire X.


Quoiqu'il en soit, il faut que tu te poses deux question:
- Qu'est ce que ça change si ta matrice n'est pas inversible ?
- Qu'est ce que ça change si ta matrice est égale à une partie fondamentale (grandes valeurs propres) non inversible + une matrice scalaire de plein rang. Cette partie scalaire est petite et provient d'un bruit, donc la matrice complète est numériquement inversible, mais sa partie fondamentale ne l'est pas

[tachmou]

Je ne vais pas m'attarder sur la théorie

Désolé d'être lourd... mais étant donné que je suis du domaine, ça m'intéresse et je pense que si tu détailles un peu plus, je pourrais vraiment t'aider. Je travaille sur des problématiques vraiment similaires, en utilisant des matrices particulièrement mal foutues, en provenance de modèles d'hélico (ce qui implique que j'utilise des matrices vraiment mal foutues - et pourtant j'ai accès à des modèles très détaillés et réalistes).

Je suis même intéressé que tu exposes le système sur lequel tu travailles, l'ordre de grandeur de tes matrices ABCD de ta représentation d'état, et qu'est-ce que tu souhaites piloter (quels paramètres de sortie?) et avec quoi (quels paramètres d'entrée?).

De plus :

Dit autrement, on a contraint les degrés de liberté de inv(X) pour le stabiliser. Une autre approximation serait de la forme X= B'*D*B + u*Identité, ou D n'est PAS de plein rang.

(...)

- Qu'est ce que ça change si ta matrice est égale à une partie fondamentale (grandes valeurs propres) non inversible + une matrice scalaire de plein rang. Cette partie scalaire est petite et provient d'un bruit, donc la matrice complète est numériquement inversible, mais sa partie fondamentale ne l'est pas

Je trouve que la première remarque que je cite de VV33D est vraiment très pertinente. La seconde peut peut-être aider aussi je pense.

Dernière remarque : la mise en place d'une loi de commande par modèle de référence pourrait peut-être t'aider... je n'ai pas trop le temps de t'expliquer comment ça fonctionne, mais en fouillant un peu sur le net tu trouveras des infos à ce sujet.

[Matersss]

Bonjour,

Le principe de la pseudo diagonalisation est expliqué un peu en pièce jointe et dans cette thèse (p323 du pdf)
[/URL]

En gros, voici mes matrices A1 et A2 (mais j'ai des matrices différentes quand je considère une autre bande de fréquences).
Le but est de trouver une matrice de découplage de degré 2 où chaque terme est sous la forme "nu1(1)+nu1(3)*s+nu1(5)*s^2" avec s=jw.

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
 
A1 =
 
  1.0e+008 *
 
   0.000037053377578   0.000013653106800   0.000056502523191  -0.000131862241888  -0.001063494798712  -0.007298866753783
   0.000013653106800   0.000013085101446   0.000068594251531   0.000087007543541   0.000093157563422  -0.003370114919682
   0.000056502523191   0.000068594251531   0.000369538836515   0.000603213592036   0.001255235431976  -0.015167571955440
  -0.000131862241888   0.000087007543541   0.000603213592036   0.002752009163470   0.011950092162123   0.014515066162804
  -0.001063494798712   0.000093157563422   0.001255235431976   0.011950092162123   0.059731915375995   0.168499138818109
  -0.007298866753783  -0.003370114919682  -0.015167571955440   0.014515066162804   0.168499138818109   1.495276029807467

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
 
A2 =
 
  1.0e+007 *
 
   0.000992742524945  -0.000725136577246   0.001501727210123  -0.002882554732437  -0.029284734754214  -0.023228797178576
  -0.000725136577246   0.000529735336306  -0.001073921249277   0.002096942596351   0.021620685702015   0.016493353650265
   0.001501727210123  -0.001073921249277   0.010022240385932  -0.007253517479930   0.033220571733667  -0.194833388487621
  -0.002882554732437   0.002096942596351  -0.007253517479930   0.009449760787805   0.056096111780362   0.127057271712980
  -0.029284734754214   0.021620685702015   0.033220571733667   0.056096111780362   1.639203376997639  -0.912018177516494
  -0.023228797178576   0.016493353650265  -0.194833388487621   0.127057271712980  -0.912018177516494   3.833926533865465

Les matrices sont carrés, symétriques et normalement définies positives par construction.
Juste en multipliant les matrices par leurs inverses respectifs, on voit qu'il y a un problème.
Pour plus de clarté, j'ai simplifié les conditions sur l'inversibilité de A1 et A2.
Je récupère le vecteur propre correspondant à la valeur propre maximale pour chacune des colonnes et ça me donne la matrice de découplage.

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
 
[V1,D1]=eig(inv(A2)*A1);
ind1=find(diag(D1)==max(diag(D1)));
nu1=V1(:,ind1);
 
[V2,D2]=eig(inv(A1)*A2);
ind2=find(diag(D2)==max(diag(D2)));
nu2=V2(:,ind2);
 
k1=[nu1(1)+nu1(3)*s+nu1(5)*s^2;nu1(2)+nu1(4)*s+nu1(6)*s^2];
k2=[nu2(1)+nu2(3)*s+nu2(5)*s^2;nu2(2)+nu2(4)*s+nu2(6)*s^2];
K=[k1,k2];

Dans certains cas, il me renvoie des vecteurs propres complexes où me renvoie des vecteurs propres qui ne découplent absolument pas mon système (et ça ne fonctionne pas en général quand un warning indique que les matrices sont mal conditionnées)

Matersss

[tachmou]

Ok, je commence à entrevoir tes soucis.

Quelques questions :
- c'est quoi A1 et A2 ? ça sort d'où ces matrices ? Comment-se fait-il qu'elles soient symétriques ? (ça, c'est particulièrement merdique d'ailleurs... sauf si les éléments principaux sont sur la diagonale... enfin tout dépend de ce que tu veux faire : d'où les questions suivantes)
- Tu utilises des représentations d'état ? Si oui, c'est quoi tes états ? Tu souhaites piloter quoi avec quoi (sorties et entrées) ?

Je n'ai pas le temps de regarder en détail la thèse que tu m'as link. De plus, comme je travaille avec des contraintes spécifiques et des critères d'optimisation particuliers, aucune des méthodes classiques ne s'appliquent à mes problématiques ==> du coup je connais mal les méthodes classiques (Hinf, mu-synthèse, etc...). Par contre, j'ai dû me poser les questions de la commandabilité de systèmes fortement couplés. C'est pourquoi je pense pouvoir t'aider. Mais il me faut encore quelques détails pour être sûr que je ne t'indique pas de mauvaises pistes.

Pour l'instant, je reste sur ma position :
- utiliser les pseudo inverses (pinv sous matlab) ou conditionner ta matrice pour n'inverser que les parties qui t'intéressent (cf méthodes de VV33D). Pour faire ça de manière "intelligente", il vaut mieux bien connaître le système d'abord (d'où mes précédentes questions).
- mettre en place une loi de commande par modèle de référence... après reste à voir s'il y a des documents accessibles à ce sujet en fait (ce sont mes encadrants qui m'ont indiqué ces méthodes) : c'est un truc que je suis en train d’approfondir pour ma thèse... mais je pense que tu dois pouvoir trouver des explications du principe à droite à gauche...

Et tout ça... c'est sans prendre en compte tes saturations actionneurs... bref : tu ne pourras pas faire des miracles si on ne te fournit pas les moyens nécessaires, bien sûr.

Cordialement,

[Matersss]

Mes matrices A1 et A2 sont construites de la façon suivante :
g1 et g2 correspondent aux réponses fréquentielles de ma matrice de transfert (je travaille avec ce type de représentation et u1->y1,u2->y2).
W sert à pondérer certaines gammes de fréquences (mais il est à 1 pour le moment).
Si j'utilise des gammes de fréquences trop espacées, les termes G(1,1) et G(1,1)*s² n'ont pas du tout la même magnitude à basse fréquence et à haute fréquence ce qui peut induire des problèmes. Mais j'ai déjà essayé de pondérer certains coefficients sans résultats.

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
    g1=[W*squeeze(freqresp(G(1,1),wt));W*squeeze(freqresp(G(1,2),wt));...
       W*squeeze(freqresp(G(1,1)*s,wt));W*squeeze(freqresp(G(1,2)*s,wt));...
       W*squeeze(freqresp(G(1,1)*s^2,wt));W*squeeze(freqresp(G(1,2)*s^2,wt))];
    g2=[W*squeeze(freqresp(G(2,1),wt));W*squeeze(freqresp(G(2,2),wt));...
       W*squeeze(freqresp(G(2,1)*s,wt));W*squeeze(freqresp(G(2,2)*s,wt));...
       W*squeeze(freqresp(G(2,1)*s^2,wt));W*squeeze(freqresp(G(2,2)*s^2,wt))];
 
  a1=real(g1);
  b1=imag(g1);
  A1=a1*transpose(a1)+b1*transpose(b1);
  a2=real(g2);
  b2=imag(g2);
  A2=a2*transpose(a2)+b2*transpose(b2);

Mon but est d'étudier sur ce problème la commande décentralisée. Donc je souhaite découpler le système et utiliser des correcteurs monovariables.

J'ai déjà essayé pinv, ça me donne une matrice bien différente de celle obtenue par inv mais j'ai toujours des matrices assez aberrantes (du moins dans certains cas) quand je fais pinv(A1)*A1

Et oui, il y aura pas mal de complications quand je vais repasser au vrai modèle

Je vais continuer sur les pistes données.
En tout cas, merci pour votre aide

Matersss

[tachmou]

Oulala, c'est compliqué tout ça...

Je ne comprends pas pourquoi tu fais tout ça...

Bon, de ce que je vois :
- tes matrices A sont symétriques, et avec des valeurs qui peuvent changer du tout au tout, sachant que les éléments principaux ne sont même pas autour de la diagonale : ça part très très mal
- tu utilises une méthode - qui est certainement très bien dans certains cas - mais je ne vois pas pourquoi tu utilises ça sachant que la suite va être très problématique
- ... alors à partir de là, j'ai encore 15000 questions...

Bref, tout ça pour dire que je pense que je vais avoir du mal à t'aider comme ça via un forum, d'autant que c'est en train de partir dans tous les sens (ce qui est un peu normal en thèse : le problème c'est que via un forum, on en a pour 15 ans pour s'en sortir... et j'ai pas trop le temps non plus).

Malgré tout, je te conseille fortement de :
1) prendre du recul (tu me sembles absorbé par la méthode que tu essayes d'appliquer : c'est super, je connais cette sensation, mais j'ai peur que tu te bloques à cause de ça, alors qu'à côté il existe peut-être des solutions bien plus appropriées)
2) jeter un coup d'oeil aux lois de commandes par modèle de référence en feedforward (les calculs théoriques peuvent être un poil hardus selon ce que tu souhaites faire... quoique les calculs que tu m'exposes me semblent plus compliqué alors ça devrait aller). J'insiste lourdement ! Apparemment, vu ce que tu souhaites faire, pour moi, c'est clair : ça répond à tous tes problèmes (découplage, contrôle fréquentiel, contrôle axe par axe, conservation d'un feedback à structure contrainte). En plus, tu peux maîtriser tes performances et la puissance de tes actionneurs, selon l'ordre de la loi que tu mets... Seule véritable contrainte : posséder un modèle réaliste de ton système de base... et ça m'a l'air d'être le cas !
3) je reste sur ma position de l'utilisation de pinv OU de conditionner ta matrice (en la bidouillant en gros) pour obtenir une matrice plus intéressante à tes yeux (ne pas hésiter à faire apparaître des zéros dans certains cas ou d'utiliser un artifice faisant apparaître des faibles valeurs dans d'autres cas) qui pourrait du coup être plus facile à inverser par exemple. Même si tu n'as pas réussi à l'utiliser efficacement dans certains cas, cela ne veut pas dire que c'est inutile : il faut juste le faire convenablement... et ce n'est pas tout le temps évident.

Le point 3) est d'ailleurs "automatiquement" utilisé dans 2) à moment donné... sauf si tu utilises une inversion dynamique auquel cas tu as encore moins de problème (et il se peut que ça corresponde d'ailleurs à ce que tu souhaites faire !!!)

Si l'occasion se présente, on peut toujours essayer de se rencontrer : ce sera beaucoup plus simple je pense...

Dans tous les cas : ce n'est pas un problème matlab que tu as, mais un problème avec ce que tu souhaites faire :-)

A la base, ici, pas grand monde peut t'aider malheureusement. VV33D est le plus à même de t'indiquer les méthodes mathématiques, et il l'a très bien fait. Et perso, je ne peux pas t'aider beaucoup plus, surtout via un forum, désolé.

D'ailleurs, j'y pense : essayes de contacter les personnes qui ont déjà appliqué la méthode que tu essayes d'appliquer... ils auront certainement déjà rencontré ces problèmes et pourront peut-être t'indiquer des pistes de solutions qui te conviendront mieux ?

[Matersss]

Je vais continuer sur les pistes données et j'ai d'autres pistes à explorer.
A la base, cela fonctionne à peu près par optimisation mais c'est juste que cette méthode permet de s'affranchir de l'optimisation et de trouver la solution optimale théoriquement. Mais c'est vrai qu'en pratique, ça m'a l'air bien compliqué.

Merci pour ces conseils Matlab, mathématiques et d'ordre général.
Matersss

[tachmou]

Salut,

Je remonte ce sujet car j'ai trouvé un excellent document qui peut peut-être aider (autour de la page 20).

A+

[Matersss]

Merci tachmou pour ce document.
Finalement, les idées proposées par tachmou ont fonctionné à moitié donc je ne vais pas poursuivre mais j'explique quand même ce que j'ai fait si ça peut aider quelqu'un :

J'ai essayé de chercher l'inverse de mes matrices en essayant de résoudre AX=B.
En posant B l'identité, cela revient à chercher X l'inverse de A.
Comme B et X sont des vecteurs, on peut décomposer le problème n fois (n étant la taille de la matrice carrée).
En gros, pour une 4*4, on résout :
AX1=[1,0,0,0]
AX2=[0,1,0,0]
AX3=[0,0,1,0]
AX4=[0,0,0,1]
Et on a X=[X1,X2,X3,X4]
Matlab possède de nombreuses routines pour trouver X sur les problèmes de type AX=B (bicg, bicgstab, cgs, gmres, lsqr,
minres, pcg, symmlq, qmr).
Cela peut améliorer la résolution de l'inverse mais visiblement seulement dans certains cas et j'ai fait un mixe des méthodes de résolution.

Sinon, d'autres idées que j'avais testées étaient de décomposer A avec qr, lu et svd et prendre l'inverse de chacune des matrices. Encore une fois, cela peut améliorer l'inverse dans certains cas.

Matersss