Discussion :

[svagrim]

Bonjour,

je cherche a résoudre un système d'équation du type Ax = B avec A et B des matrices.
et je voudrais spécifier que x ne peux avoir que des valeurs positives.
il y a un moyen de faire ça facilement ?

merci

[Cheetor]

Il y a la possibilité suivante :

Code :

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

1
2
opts.POSDEF = true;
X = linsolve(A,B,opts);

De manière plus globale, la page suivante liste toutes les options que tu peux ajouter : http://www.mathworks.fr/help/matlab/ref/linsolve.html

[svagrim]

ok merci.
j'avais vu cette option mais j'avais pas compris ce qu'ils voulaient dire par posdef vu que c'est dit que c'est des options en rapport avec la matrice A et non x.

j'ai essayé et j'ai ce message d'erreur :

Code :

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

1
2
Error using linsolve
the option slected by the combination of fields in the structure array is currently not supported

si j'essaye avec les autres options de linsolve (UT etc ) j'ai le meme message d'erreur.

à quoi c'est du ?

[shaiHulud]

Bonjour,

Pas sur de moi, mais a priori, posdef signifie définie positif et porte sur A (rien a voir avec la positivité des éléments de x ! - revoir la différence entre une matrice symétrique définie positive et une matrice dont les éléments sont positifs)

Une solution: passer par un optimiseur plus général comme fmincon, qui permet de passer des contraintes d'égalité, d'inégalité et/ou d'inégalité coordonnée à coordonnée sur la variable optimisée x.

je suis surpris que ces contraintes ne puissent pas être passées à linsolve, mais cela doit être du aux factorisations utilisées que l'on ne peut pas naturelement contraindre de la sorte

[svagrim]

du coup avec fmincon, si je veux résoudre par example:
Je cherche C = [C1, C2, C3};
je mesure:
X = [X1, X2, X3];
et je connais: a = [A1 A2 A3; A4 A5 A6; A7 A8 A9];

et ça doit faire:
A1*C1 + A2*C2 + A3*C3 = X1;
A4*C1 + A5*C2 + A6*C3 = X2;
A7*C1 + A8*C2 + A9*C3 = X3;
Avec C1>0, C2>0, C3>0;
C1+C2+C3 <= 1;

du coup avec fmincon je fais:

Code :

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

1
2
function f = myfunc(x,A);
f= A*x;

je met A dans la fonction parce que j'ai beaucoup de A différent et je dois les "injecté" automatiquement dans l'équation

Code :

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

1
2
3
4
 
x0 = [0.5 0.5 0.5]
 
[x,fval] = fmincon(@myfun,x0,[-1 0 0; 0 -1 0; 0 0 -1; 1 1 1];[0; 0; 0; 1]);

pour le moment je comprend comme ça fmincon mais je met ou X et A ?

[le fab]

salut

par exemple :

Code :

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

[x,fval] = fmincon(@(x) myfun(x,A),x0,[-1 0 0; 0 -1 0; 0 0 -1; 1 1 1];[0; 0; 0; 1]);

fabien

[svagrim]

merci,

j'ai essayé avec 2 inconnus au lieu de droit pour le moment, et ça me sort le message d'erreur:

Code :

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

1
2
3
4
Warning: The default trust-region-reflective algorithm does not solve problems with the constraints you have specified. FMINCON will use the active-set algorithm instead. For information on applicable algorithms, see Choosing the Algorithm in the documentation.
 
Error using fmincon (line 708)
User supplied objective function must return a scalar value.

ma ligne de code était:

Code :

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

1
2
 x0 = [0.5;0.5];
 [concentration, fval] = fmincon(@(mu_measured) linear_matrix_eq(mu_measured,coeff_attenuation),x0,[-1 0; 0 -1; 1 1],[0; 0; 1]);

avec les valeurs:

mu_measured =

   0.0218
   0.1020

coeff_attenuation =

    0.1984    0.0066
    0.3463    0.0187

je comprend pas le coup du scalar value.

[le fab]

ta fonction doit retourner un nombre (scalar) et non une matrice (un vecteur dans ton cas)
en effet, que veut dire "minimiser un vecteur" ? qui est le plus petit de [1 2] et de [3 1] ?
en general on va minimiser la norme
ex

Code :

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

[x,fval] = fmincon(@(x) norm(myfun(x,A)),...)

[svagrim]

vui c'est ça.
pour moi une resolution matricielle tu cherches le minimum pour la norme dans la plupart des cas simple.

genre t'obtiens le vecteur "T"
A1*C1 + A2*C2 + A3*C3 - X1 = T(1);
A4*C1 + A5*C2 + A6*C3 - X2 = T(2);
A7*C1 + A8*C2 + A9*C3 - X3 = T(3);

et ||T|| doit etre minimiser.

je pensais que la fonction le "ferait tout seul" vu que c'est la type de résolution de système d'équation matricielle basique.

du coup je vais essayer en ajoutant comme tu dis.

[svagrim]

en fait me manquait C aussi dans ma fonction ... je me suis mélanger les pinceaux...

du coup j'ai fait:

Code :

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

1
2
3
function fmulti = linear_matrix_eq(x,AA,C)
 
fmulti = AA*C-x;

et

Code :

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

1
2
 x0 = [0.5;0.5];
[concentration, fval] = fmincon(@(conc) norm(linear_matrix_eq(mu_measured,coeff_attenuation,conc)),x0,[-1 0; 0 -1; 1 1],[0; 0; 1]);

ça tourne. ça me ressort des chiffres qui, même si je suis pas hyper convaincu pour le moment, ne sont pas non plus totalement absurdes.

par contre à chaque résolution il me display ça:

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
Warning: The default trust-region-reflective algorithm does not solve problems with the constraints you have specified. FMINCON will use the active-set
algorithm instead. For information on applicable algorithms, see Choosing the Algorithm in the documentation.
> In fmincon at 504
  In multi_engergy_processing_guide>Validate_button_Callback at 400
  In gui_mainfcn at 96
  In multi_engergy_processing_guide at 42
  In @(hObject,eventdata)multi_engergy_processing_guide('Validate_button_Callback',hObject,eventdata,guidata(hObject))
Warning: Your current settings will run a different algorithm (interior-point) in a future release.
> In fmincon at 509
  In multi_engergy_processing_guide>Validate_button_Callback at 400
  In gui_mainfcn at 96
  In multi_engergy_processing_guide at 42
  In @(hObject,eventdata)multi_engergy_processing_guide('Validate_button_Callback',hObject,eventdata,guidata(hObject))
 
Local minimum possible. Constraints satisfied.
 
fmincon stopped because the predicted change in the objective function
is less than the default value of the function tolerance and constraints
are satisfied to within the default value of the constraint tolerance.
 
<stopping criteria details>
 
Active inequalities (to within options.TolCon = 1e-06):
  lower      upper     ineqlin   ineqnonlin

ça m'a l'air plutot positif comme message. par contre j'ai 1000x1000x1000 fois où je fais ce calcul. alors s'il pouvait ne pas affiché ça me ferait gagné énormément de temps de calcul ....
est-ce que vous sauriez comment l'empêcher d'afficher tout ça ?

[le fab]

il suffit de définir l'algorithme 'interior point' au lieu de celui par defaut pour virer le warning et quand au resultat de l'optimisation il suffit de mettre la propriété 'display' à "off"
le tout avec optimset

[svagrim]

ca a l'air de marché maintenant.

je veux deux trois tests pour vérifier mais ça a l'air d'aller.

pour contre j'ai mis [],[],[],[],[]
pour Aeq, beq, lp, up, nonlcon
c'est pas grave ?

j'ai l'impression que Aeq et beq, lp et up dans mon cas j'aurais pu les utiliser au lieu de faire minimiser la norme et de dire avec les matrices que x doit etre positif nan ?

et par contre, gros gtros problème.
j'estime le temps que ça va prendre à 900 heures à peu prêt ....
j'ai 875*875*500 fois ces equations à résoudre et ça prend juste un temps phénoménale.
il y aurait un moyen d'accélérer tout ça ?

j'ai mis aussi dans l'optimization

Code :

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

'MaxIter',600,'TolCon', 1e-3, 'TolFun', 1e-3, 'TolX', 1e-5

mais ça n'a quasiment pas changé le temps de calcul.

comme mes équations linéaires il n'y a pas un truc que je pourrais simplifier quelque part pour accélérer les calculs ? (j'ai essayé 'UseParallel', 'Always' mais ça n'a fait strictement aucune différence).

[le fab]

si, il faut utiliser au maximum les contraintes (A, b, Aeq, bEq ..)
plus tu contraint ton pb, plus la résolution sera rapide
d'autant que pour le moment tu n'a pas gerer la contrainte "ts les élément de x positifs"

[svagrim]

hum, je comprend plus rien ...

que mes x doivent etre positifs c'est définit par (ce qui est souligné):

Code :

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

[concentration, fval] = fmincon(@(conc) norm(linear_matrix_eq(mu_measured,coeff_attenuation,conc)),x0,[-1 0; 0 -1; 1 1],[0; 0; 1]);

et ça marche. ça me ressort que des résultats entre 0 et 1 conformément à ce que je souhaitais.

mon égalité est défini par ma fonction vu que je cherche a minimiser la différence "fonction - vecteur résultat" du coup je comprend pas trop ce que je pourrais ajouter dans Aeq, beq, lb, ub.

en fait si j'utilise Aeq, beq ma fonction ne sert à rien.
je suis pas sur que ce soit la bonne fonction en fait fmincon.
elle a l'air bien plus complexe que ce que j'ai besoin non ?
mon problème (message 5) n'est pas résolvable par une autre fonction ? (qui serait peut etre plus rapide en plus du coup).

[le fab]

au temps pour moi j'ai lu un peu vide
fmincon est peut être un peu complexe mais c'est la fonction d'optimisation lorsque le problème est contraint

cela dit tes contraintes très simples (bornes) devraient être traité avec lb, ub et non A, b
ces contraintes simples sont également prisent en compte par fminbnd
(mais cela ne sera pas forcement plus rapide)

après si ton pb est quadratique, quadprog sera très rapide
quand à fmincon, si tu sais calculer les matrices hessian de ta fonction, ca ira beaucoup plus vite aussi

enfin tu peux egalement jouer sue les propriétés 'FinDiffRelStep' et 'TypicalX'

[svagrim]

merci pour les précisions.

le problème n'est pas quadratique.

je vais du coup mettre mes bornes comme tu dis au lieu de passer par A et b.

je vais voir pour les hessiens, j'imagine que c'est faisable oui y a pas trop de raison que j'y arrive pas, ma fonction est assez basique. faut juste que j'aille rechoper la définition exacte quelque part.

sinon j'ai vu qu'il y avait la fonction lsqlin qui pourrait correspondre aussi sur matlab.
je pourrais peut etre tenter ça aussi pour voir.

je vous tiens au courant.
en tout cas merci pour l'aide !

[le fab]

je vais voir pour les hessiens, j'imagine que c'est faisable oui y a pas trop de raison que j'y arrive pas, ma fonction est assez basique. faut juste que j'aille rechoper la définition exacte quelque part.

si tu n'y arrive pas, j'ai la symbolic math toolbox, je peux te filer un coup de main

[svagrim]

merci pour la proposition !

j'ai pas pu travailler dessus aujourd'hui mais je m'y met demain.
si j'y arrive pas je te tiens au courant.

[svagrim]

je vais peut etre avoir l'air bete la mais j'ai regardé un hessian et bon c'est juste des derivées partielles. rien qui fasse particulièrement peur.

sauf que moi ma fonction est hyper basique.
si je dérive une vois j'ai des constantes, je dérive une seconde fois j'ai des zéros partout ...

"Je cherche C = [C1, C2, C3];
je mesure:
X = [X1, X2, X3];
et je connais: a = [A1 A2 A3; A4 A5 A6; A7 A8 A9];

et ça doit faire:
A1*C1 + A2*C2 + A3*C3 = X1;
A4*C1 + A5*C2 + A6*C3 = X2;
A7*C1 + A8*C2 + A9*C3 = X3;
Avec C1>0, C2>0, C3>0;
C1+C2+C3 <= 1;"

du coup je fais quoi je précise une matrice [0 0 0; 0 0 0 ; 0 0 0] ?

[le fab]

je pennse oui
essaye et regarde ce que ca donne