Discussion :

[cpalperou]

Bonjour,
voici mon problème: j'ai 2 vecteurs: un vecteur temps timeS et un vecteur représentant l'altitude pour chacun de ces instants altS (tout 2 de dimensions identiques = 1*80)
Je dois faire un fit de l'altitude en fonction du temps en ne prenant que les 20 premières valeurs. Je fais donc:

Code :

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

1
2
 
altFIT=polyval(polyfit(timeS(1:20),altS(1:20),2),altS)

Cela fonctionne très bien mais, comme c'est le résultat d'un moindre carré, altFIT est parfois inférieur, parfois supérieur à altS. Cependant, je dois contraindre mon problème: altFIT doit être toujours inférieur à altS.
Comment faire ceci?
Merci

[shaiHulud]

Bonjour,

Je ne pense pas que polyfit le permette.
Regarde du coté de fmincon, c'est la fonction matlab la plus générique pour faire de l'optimisation sous contrainte.

[cpalperou]

Merci pour votre réponse, je vais essayer de comprendre cette fonction (c'est pas gagné!!!)

[shaiHulud]

En fait, elle permet de contraindre facilement les paramètres. Mais toi tu cherches à contraindre les valeurs, il faut donc utiliser une contrainte non linéaire d'inégalité.

Quel est l’intérêt d'avoir un fit inférieur aux données ??

[cpalperou]

Eh bien, je n'y arrive pas!!!
Pour bien comprendre, dans un premier temps, j'essaie sans contrainte. J'utilise donc fminunc en utilisant comme condition initiale, le résultat de polyfit (ce n'est peut etre pas une bonne idée?)

Code :

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

1
2
3
4
5
6
7
8
9
 
function H=FitAlt(alt,time)
a0=polyfit(time,alt,2);
polH=fminunc(@fitness,a0,[],[])
H=polH(3)+polH(2)*time+polH(1)*time^2
%-------------------------------------------------------
function y=fitness(altSP,t,a)
Hfit=a(3)+a(2)*t+a(1)*t.^2;
y=altSP-Hfit

et dans la fenêtre de commande:

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
L=find(altSP(:,1)>100);
alt=altSP(L,1)';time=timeSP(L);
H=FitAlt(alt,time)
Error using FitAlt>fitness (line 19)
Not enough input arguments.
 
Error in fminunc (line 254)
        f = feval(funfcn{3},x,varargin{:});
 
Error in FitAlt (line 12)
polH=fminunc(@fitness,a0,[],[])
 
Caused by:
    Failure in initial user-supplied objective function evaluation. FMINUNC cannot continue.

[shaiHulud]

* Il faut mettre fitness au carré (minimiser l'écart au carré)
* La fonction que tu passes à fmincon doit avoir un seul argument (celui sur lequel on minimise). Il faut lui passer @(a) fitness(altSP,t,a)

[cpalperou]

En fait, elle permet de contraindre facilement les paramètres. Mais toi tu cherches à contraindre les valeurs, il faut donc utiliser une contrainte non linéaire d'inégalité.

Quel est l’intérêt d'avoir un fit inférieur aux données ??

Ouh là!!! comprends pas!!

L'intérêt répond à une contrainte physique (en rapport avec la réfraction)

[cpalperou]

Mille merci pour votre aide, mais ... je n'y arrive toujours pas. J'ai un epu modifié mon code:

Code :

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

1
2
3
4
5
6
7
8
9
10
11
function H=FitAlt(alt,time)
a0=polyfit(time,alt,2);
polH=fminunc(@(a) fitness(alt,time,a),a0,[],[])
H=polH(3)+polH(2)*time+polH(1)*time^2
%--------------------------------
function y=fitness(altSP,t,a)
%     t=time;
%     altSP=alt;
    Hfit = a(3)+a(2)*t+a(1)*t.^2;
 
    y=(altSP-Hfit)^2;

et j'ai cette erreur:

Code :

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

1
2
3
4
5
6
7
8
9
10
11
Error using FitAlt>@(a)fitness(alt,time,a)
Too many input arguments.
 
Error in fminunc (line 254)
        f = feval(funfcn{3},x,varargin{:});
 
Error in FitAlt (line 12)
polH=fminunc(@(a) fitness(alt,time,a),a0,[],[])
 
Caused by:
    Failure in initial user-supplied objective function evaluation. FMINUNC cannot continue.

[shaiHulud]

Code :

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

L'intérêt répond à une contrainte physique (en rapport avec la réfraction)

Je veux bien un peu plus de précision!

Code :

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

Ouh là!!! comprends pas!!

En gros, avec la fonction fitness que tu as posté, il faut minimiser @(a) fitness(altSP,t,a).^2, sous contrainte que fitness>=0.

Cette contrainte sera injectée via le 9eme argument <nonlcon> de fmincon. Cette fonction nonlcon doit renvoyer deux outputs : 0 et -fitness. 0 correspond à une égalité (tu n'en as pas, donc 0=0), et le 2eme output à une contrainte<=0

[cpalperou]

Je veux bien un peu plus de précision!

En voici:
soit un instrument satellitaire en orbite autour de la terre. Cet instrument vise une étoile à travers l'atmosphère. La ligne de visée est déviée vers le bas du fait de la réfraction atmosphérique. L'instrument nous fourni l'altitude de chaque ligne de visée en fonction du temps. L'altitude d'une ligne de visée étant défini comme la distance minimale entre la ligne de visée et la surface terrestre.
Si on enlève l'atmosphère, la ligne de visée n'est plus déviée et arrive en ligne droite avec une nouvelle altitude. Cette altitude est nécessairement inférieur à celle avec atmosphère. Je l'obtiens en faisant un fit de l'altitude (avec réfraction) fourni par l'instrument.

[shaiHulud]

Merci pour les détails. Dans ton cas, j'ai l'impression qu'une approche par quantile serait préférable à une approche par moindre carrés.
As-tu des informations sur la loi de tes perturbations ? Car un MCO contraint, ou des quantiles, vont surement sur estimer ta vraie courbe en saturant la contrainte, alors que la vraie courbe ne la sature surement jamais.

[cpalperou]

Voici ce que je cherche à faire. En rouge, c'est l'altitude en fonction du temps, celle que je dois fitter. Et en bleu, le résultat de polyval(polyfit(time(L),alt(L),2),time); L servant à sélectionner les altitudes > 100 km

Il ne faut pas que le fit soit supérieur. C'est tout simple et pourtant, je n'arrive pas à le programmer!

[cpalperou]

Cette contrainte sera injectée via le 9eme argument <nonlcon> de fmincon. Cette fonction nonlcon doit renvoyer deux outputs : 0 et -fitness. 0 correspond à une égalité (tu n'en as pas, donc 0=0), et le 2eme output à une contrainte<=0

donc, je crée une fonction nonlcon. Mais avec quel paramètre d'entrée? Je suis paumé!

[shaiHulud]

je n'arrive pas à le programmer!

fmincon devrait marcher. Ou bloques-tu ?

[shaiHulud]

nonlcon, tout comme ta fonction fitness, prends en argument tout le monde (les data et les paramètres).
Par contre, fmincon s'attends à avoir une fonction ne dépendant que des paramètres (tout comme pour sa fonction objectif).
Il faudra donc passer @(a) nonlcon(altSP,t,a) à fmincon

[cpalperou]

Merci encore de ton aide! Je me sens nul, j'ai un peu honte là!! J'ai du mal avec ces fonctions matlab qui font intervenir des @ !!
Voici ou j'en suis: j'ai une première fonction FitAlt.m:

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
13
 
function H=FitAlt(alt,time)
 
a0=polyfit(time,alt,2);
polH=fmincon(@(a)fitness(alt,time,a).^2,a0,[],[],[],[],[],[],@(a)nonlcon(alt,time,a));
 
H=polH(3)+polH(2)*t+polH(3)*t.^2;
 
%--------------------------------
function y=fitness(altSP,t,a)
 
    Hfit = a(3)+a(2)*t+a(1)*t.^2;
    y=(altSP-Hfit).^2;

et une deuxième fonction pour la contrainte: nonlcon.m

Code :

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

1
2
3
4
function [c,ceq]=nonlcon(altSP,t,a)
 
c=0;
ceq=-fitness(altSP,t,a)

puis je lance le programme:

Code :

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

H=FitAlt(alt,time)

et j'ai l'erreur:

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
13
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 FitAlt at 8 
Warning: Your current settings will run a different algorithm (interior-point) in a future release. 
> In fmincon at 509
  In FitAlt at 8 
Error using fmincon (line 708)
User supplied objective function must return a scalar value.
 
Error in FitAlt (line 8)
polH=fmincon(@(a)fitness(alt,time,a).^2,a0,[],[],[],[],[],[],@(a)nonlcon(alt,time,a));

[shaiHulud]

Les @ permettent d'utiliser les fonctions partielles (fonctions de a, avec les autres arguments ayant des valeurs fixés= à tes données)

Error using fmincon (line 708)
User supplied objective function must return a scalar value.

C'est plutot clair, ta fonction objectif (celle qu'on minimise) doit renvoyer un réel, pas un vecteur.
Ici, on a oublié de sommer les erreurs ! Il faut bel et bien minimiser la somme des erreurs aux carré.

polH=fmincon(@(a) sum(fitness(alt,time,a).^2),a0,[],[],[],[],[],[],@(a)nonlcon(alt,time,a));

Quand au warning, tu peux passer optimset('Algorithm','interior-point') en dernier argument de fmincon pour choisir le bon algo des le début.

Autre point important, ton point initial a0 n'est pas une idée stupide, mais il ne vérifie pas les contraintes, donc l'optimiseur va raller ! Il faut que tu en trouves un autre.

[cpalperou]

Après correction de mon programme (j'ai inclus la somme des carrés dans la fonction fitness), toujours des problèmes:
le programme:

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
 
function H=FitAlt(alt,time)
a0=polyfit(time,alt,2)
a0(3)=a0(3)-5; % pour etre sur
polH=fmincon(@(a)fitness(alt,time,a),a0,[],[],[],[],[],[],@(a)nonlcon(alt,time,a));
 
H=polH(3)+polH(2)*time+polH(1)*time.^2;
 
%--------------------------------
function y=fitness(altSP,t,a)
 
    Hfit = a(3)+a(2)*t+a(1)*t.^2;
 
    y=sum((altSP-Hfit)).^2;

Je n'ai pas changé lnoncon.m
et j'ai encore une erreur:

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 FitAlt at 9 
Warning: Your current settings will run a different algorithm (interior-point) in a future release. 
> In fmincon at 509
  In FitAlt at 9 
Undefined function 'fitness' for input arguments of type 'double'.
 
Error in nonlcon (line 4)
ceq=-fitness(altSP,t,a)
Error in FitAlt>@(a)nonlcon(alt,time,a) (line 9)
polH=fmincon(@(a)fitness(alt,time,a),a0,[],[],[],[],[],[],@(a)nonlcon(alt,time,a));
 
Error in fmincon (line 722)
        [ctmp,ceqtmp] = feval(confcn{3},X,varargin{:});
 
Error in FitAlt (line 9)
polH=fmincon(@(a)fitness(alt,time,a),a0,[],[],[],[],[],[],@(a)nonlcon(alt,time,a));
 
Caused by:
    Failure in initial user-supplied nonlinear constraint function evaluation. FMINCON cannot continue.

[shaiHulud]

j'ai inclus la somme des carrés dans la fonction fitness

il ne faut pas, il faut l'inclure dans le @ passé à fmincon. La raison est que ta nonlcon utilise fitness (et a besoin des chaque terme, sans carré ni somme)

Il faut que tu définisses une fonction nonlcon. Si tu veux que nonlcon utilise la fonction fitness, il faut la mettre dans le même fichier.

[cpalperou]

Hallelujah !! ... ou presque
Je n'arrive pas à imposer la contrainte (j'ai des résidus positifs et d'autres négatifs). Voici mes 2 fonctions:

Code :

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

1
2
3
4
5
6
7
8
9
10
function [H,polH]=FitAlt(alt,time)
a0=polyfit(time,alt,2)
a0(3)=a0(3)-1;
options = optimset('algorithm','interior-point')
polH=fmincon(@(a)sum(fitness(alt,time,a)).^2,a0,[],[],[],[],[],[],@(a)nonlcon(alt,time,a),options);
H=polH(3)+polH(2)*time+polH(1)*time.^2;
%--------------------------------
function y=fitness(altSP,t,a)
    Hfit = a(3)+a(2)*t+a(1)*t.^2;
    y=altSP-Hfit;

et:

Code :

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

1
2
3
4
5
6
function [c,ceq]=nonlcon(altSP,t,a)
c=0;
ceq=-fitness(altSP,t,a);
function y = fitness(altSP,t,a)
Hfit = a(3)+a(2)*t+a(1)*t.^2;
y=altSP - Hfit;

Et dans la fenêtre de commande matlab, je tape [H,polH]=FitAlt(alt,time);

là, matlab m'affiche la liste des options et plusieurs fois ce même message d'avertissement:

Code :

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

1
2
3
4
5
6
7
8
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =  1.563323e-22. 
> In /opt/matlab/toolbox/optim/optim/private/backsolveSys.p>backsolveSys at 18
  In /opt/matlab/toolbox/optim/optim/private/solveAugSystem.p>solveAugSystem at 23
  In /opt/matlab/toolbox/optim/optim/private/leastSquaresLagrangeMults.p>leastSquaresLagrangeMults at 28
  In /opt/matlab/toolbox/optim/optim/private/nlpStopTest.p>nlpStopTest at 125
  In /opt/matlab/toolbox/optim/optim/barrier.p>barrier at 697
  In fmincon at 900
  In FitAlt at 11

puis ce message

Code :

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

1
2
3
4
5
6
7
No feasible solution found.
 
fmincon stopped because the size of the current step is less than
the default value of the step size tolerance but constraints are not
satisfied to within the default value of the constraint tolerance.
 
<stopping criteria details>