Discussion :

[Eleasias]

bonjour,
Je cherche à fitter un modèle sur un jeu de données. Mais ça ne marche pas bien...
J'ai écrit un code qui définit mon modèle (l'une ou l'autre des equadiff MyModel1 ou MyModel2), calcule un y_calc, fait la somme des carrés des différences avec mes valeurs mesurées (Val) et utilise cette valeur pour ajuster mes deux paramètres : k=[u, Ks]

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
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
clc
clear
funcprot(0)
 
y_exp(:,1) = [8000 23000 61000 101000 116000 192000 183000 196000]';
y_exp(:,2) = [13000 9000 53000 77000 116000 205000 232000 261000]';
y_exp(:,3) = [6000 26000 46000 67000 136000 192000 236000 303000]';
y_exp(:,4) = [8000 19000 8000 11000 8000 9000 17000 20000]';
y_exp(:,5) = [7000 16000 8000 20000 33000 40000 32000 85000]';
y_exp(:,6) = [9000 16000 15000 16000 36000 74000 57000 79000]';
        // data from 12-03-13
y_exp(:,7) = [11000 9000 51000 104000 128000 182000 192000 284000]';
y_exp(:,8) = [9000 7000 14000 26000 46000 100000 114000 162000]';
y_exp(:,9) = [11000 13000 74000 185000 226000 356000 228000 428000]';
y_exp(:,10) = [11000 24000 22000 30000 78000 112000 168000 156000]';
y_exp(:,11) = [13000 27000 111000 201000 222000 272000 316000 332000]';
y_exp(:,12) = [14000 25000 150000 260000 254000 364000 396000 452000]';
        // data from 27-03-13
y_exp(:,13) = [10000 10000 16000 33000 46000 84000 75000 104000]';
y_exp(:,14) = [7000 8000 40000 66000 135000 228000 180000 204000]';
y_exp(:,15) = [8000 12000 38000 78000 129000 180000 162000 248000]';
y_exp(:,16) = [6000 10000 11000 14000 11000 28000 23000 27000]';
y_exp(:,17) = [6000 7000 9000 25000 23000 31000 50000 49000]';
y_exp(:,18) = [5000 7000 22000 32000 42000 50000 55000 65000]';
        // data from 11-04-13
y_exp(:,19) = [7000 24000 116000 188000 224000 332000 328000 392857]';
y_exp(:,20) = [8000 15000 39000 118000 140000 196000 288000 427429]';
y_exp(:,21) = [8000 12000 43000 103000 144000 192000 300000 480857]';
y_exp(:,22) = [7000 9000 30000 39000 82000 128000 160000 232571]';
y_exp(:,23) = [5000 9000 25000 51000 75000 136000 196000 238857]';
y_exp(:,24) = [8000 11000 27000 45000 99000 148000 184000 257714]';
 
Fluxi = [148.27 148.71 140.59 51.65 54.36 52.76 148.27 148.71 140.59 233.60 238.9 244.0 148.27 148.71 140.59 51.65 47.12 43.29 241.30 241.80 242.50 ...
141.93 160.74 159.9]';
 
DOi = [0.425 0.416 0.414 1.286 0.374 0.191 1.325 0.399 0.149 1.358 0.363 0.221 0.384 0.893 0.379 0.194 0.917 0.382 0.855 0.25 0.892 0.885 0.87 0.887]';
 
//Jb
function [v] = Jb(t, y, k, j, DOi, Fluxi)
    g=12.789;
    h=0.000196202;
    i=54.7014;
    v = g + h*y+ i*DOi(j);
endfunction
 
//Flux
function [z] = Flux(t, y, k, j, DOi, Fluxi)
    a=0.902623;
    b=0.02971;
    c=0.0005421;
    d=0.00000557133;
    e=0.0000000299771;
    f=0.0000000000653467;
    v=Jb(t, y, k, j, DOi, Fluxi);
    z=Fluxi(j)*(a-b.*v+c.*v.^2-d.*v.^3+e.*v.^4-f.*v.^5);
endfunction
 
// Model jour
function [ydot] = mymodel1(t, y, k, j, DOi, Fluxi)
        z = Flux(t, y, k, j, DOi, Fluxi);
        ydot = k(1).*y.*(z./(z+k(2)));
endfunction
 
//Model nuit
function [ydot] = mymodel2(t, y, k, j, DOi, Fluxi)
        ydot=0;
endfunction
y_calc=[];
//Least squares
function [val] = L_Squares(k)
val=[];
u = k(1);
Ks = k(2);
for j=1:24;
    y_temp=[]
    t = 0;
    while t<14;
        temp = (t-int(t));
        if (temp<0.25) then
            if int(t)==0 then
                y0 = y_exp(1,j);
                //printf('y0-i -- %g\n',y0)
                ta=[0:pas:0.25];
                y=ode(y0,0,ta,mymodel1)
            else
                y0 = y_trans3(1,$);
                //printf('y0-1 -- %g\n',y0)
                ta=[pas:pas:0.25];
                y=ode(y0,0,ta,mymodel1)
            end
            y_trans1=y;
            t=t+0.25;
        elseif (temp>=0.25)&(temp<0.75) then
            tb=[pas:pas:0.5];
            y0 = y_trans1(1,$);
            //printf('y0-1 -- %g\n',y0)
            y=ode(y0,0,tb,mymodel2)
            y_trans2=y;
            t=t+0.5
        elseif (temp>=0.75) then
            tc=[pas:pas:0.25];
            y0 = y_trans2(1,$);
            //printf('y0-1 -- %g\n',y0)
            y=ode(y0,0,tc,mymodel1)
            y_trans3=y;
            y_trans=[y_trans1'; y_trans2'; y_trans3'];
            t=t+0.25
            y_temp=[y_temp;y_trans];
        end
    end
    y_calc(:,j)=y_temp;
    for m=1:8;
        n=(m*2/pas+1-2/pas);
        y_transb(m,j)=y_calc(n,j);
    end
    resid = (y_transb(:,j) - y_exp(:,j)).*(y_transb(:,j) - y_exp(:,j)); //matrice of residues
    resid = sum(sum(resid));//residue for the experiment j
    Residus(1,j)=resid(1,1);
end
    val = sum(sum(Residus));//global residue
//    printf('val -- %g\n',val)
//    val = 1e12
endfunction
 
// starting guess
u = 0.25;
Ks = 50;
k = [u;Ks];
pas=0.01
//---------------------------------------------------------------------------
tic
options = optimset("TolFun",1D-04,"TolX",1.D-04,"MaxFunEvals",1.0D+06,"MaxIter",1.0D+06);
[param,L_Squares,flag,details] = fminsearch(L_Squares,k,options);
elapsed_time = toc();
parametres = param';
u = parametres(1); Ks = parametres(2);
//clc
parametres
if flag == 1 then
printf('Le système converge\n\n');
else
printf('Le système ne converge pas\n\n');
end;
printf('nombre d itération ---------------- %g\n',details.iterations);
printf('Somme des carrés ----------------- %g \n',L_Squares);
printf('Nombre d appels du modèle ------- %g\n',details.funcCount);
printf('Algorithme de minimisation------------ %s\n',details.algorithm);
printf('µ ');
write(%io(2),u);
printf('Ks :');
write(%io(2),Ks);
if elapsed_time<60 then
printf('Temps écoulé -------------------- %g s \n\n',elapsed_time);
else
printf('Temps écoulé -------------------- %g min \n\n',elapsed_time);
end;

Lorsque je lance ce code, la console me renvoie un "index invalid" que je ne comprend pas, précédé de quelques lignes (voir ci-après)

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
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
 lsoda--  caution... t (=r1) and h (=r2) are
     such that t + h = t at next step
      (h = pas). integration continues
      where r1 is :   0.1910482231662D+00   and r2 :   0.9882870128075D-17      
lsoda--  caution... t (=r1) and h (=r2) are
     such that t + h = t at next step
      (h = pas). integration continues
      where r1 is :   0.1910482231671D+00   and r2 :   0.1364457904114D-16      
lsoda--  caution... t (=r1) and h (=r2) are
     such that t + h = t at next step
      (h = pas). integration continues
      where r1 is :   0.1910482231671D+00   and r2 :   0.1364457904114D-16      
lsoda--  caution... t (=r1) and h (=r2) are
     such that t + h = t at next step
      (h = pas). integration continues
      where r1 is :   0.1910482231671D+00   and r2 :   0.1364457904114D-16      
lsoda--  caution... t (=r1) and h (=r2) are
     such that t + h = t at next step
      (h = pas). integration continues
      where r1 is :   0.1910482231671D+00   and r2 :   0.2102081432332D-17      
lsoda--  caution... t (=r1) and h (=r2) are
     such that t + h = t at next step
      (h = pas). integration continues
      where r1 is :   0.1910482231671D+00   and r2 :   0.2102081432332D-17      
lsoda--  caution... t (=r1) and h (=r2) are
     such that t + h = t at next step
      (h = pas). integration continues
      where r1 is :   0.1910482231671D+00   and r2 :   0.2102081432332D-17      
lsoda--  caution... t (=r1) and h (=r2) are
     such that t + h = t at next step
      (h = pas). integration continues
      where r1 is :   0.1910482231671D+00   and r2 :   0.1023537392927D-16      
lsoda--  caution... t (=r1) and h (=r2) are
     such that t + h = t at next step
      (h = pas). integration continues
      where r1 is :   0.1910482231671D+00   and r2 :   0.1023537392927D-16      
lsoda--  at t (=r1), mxstep (=i1) steps   
needed before reaching tout
      where i1 is :        500                                                  
      where r1 is :   0.1910482231711D+00                                       
Attention: Le résultat est peut être inexact.
 
lsoda--  at t (=r1), mxstep (=i1) steps   
needed before reaching tout
      where i1 is :        500                                                  
      where r1 is :   0.1048191146817D-02                                       
Attention: Le résultat est peut être inexact.
 
 !--error 21 
Index invalide.
 
at line      45 of function  called by :  
at line       2 of function  called by :  
at line      59 of function optimbase_function called by :  
at line     122 of function neldermead_variable called by :  
at line       6 of function neldermead_algo called by :  
at line      15 of function neldermead_search called by :  
at line      90 of function fminsearch called by :  
[param,L_Squares,flag,details] = fminsearch(L_Squares,k,options);
at line     133 of exec file called by :    
exec('\\Re-data\data\gere\1-Projets\ALGOVALO\Manips\Modélisation\LeastSquares_Manip_Couleur-Lumière_Monod_12-12_4.sci', -1)

J'ai fait affiché la valeur de Val et elle est toujours égale à 1 seule nombre et je n'ai pas la main sur le k (seulement l'initialisation).
Est-ce que quelqu'un y comprend quelque chose?
Merci beaucoup de l'aide que vous pourriez m'apporter.
Cordialement,
Eleasias

[Eleasias]

j'ai essayé de changer mes valeurs de u et de Ks (0.2 et 70) et la réponse est la suivante :

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
25
26
27
28
29
30
31
32
33
34
35
36
37
38
 lsoda--  caution... t (=r1) and h (=r2) are
     such that t + h = t at next step
      (h = pas). integration continues
      where r1 is :   0.6626113833734D-01   and r2 :   0.5016440108464D-17      
lsoda--  caution... t (=r1) and h (=r2) are
     such that t + h = t at next step
      (h = pas). integration continues
      where r1 is :   0.6626113833734D-01   and r2 :   0.5016440108464D-17      
lsoda--  caution... t (=r1) and h (=r2) are
     such that t + h = t at next step
      (h = pas). integration continues
      where r1 is :   0.6626113833734D-01   and r2 :   0.5016440108464D-17      
lsoda--  at t (=r1), mxstep (=i1) steps   
needed before reaching tout
      where i1 is :        500                                                  
      where r1 is :   0.6626113834290D-01                                       
Attention: Le résultat est peut être inexact.
 
lsoda--  at t (=r1), mxstep (=i1) steps   
needed before reaching tout
      where i1 is :        500                                                  
      where r1 is :   0.6261155299418D-02                                       
Attention: Le résultat est peut être inexact.
 
lsoda--  neq (=i1) .lt. 1     
      where i1 is :          0                                                  
 !--error 9999 
illegal input
at line      28 of function  called by :  
at line       2 of function  called by :  
at line      59 of function optimbase_function called by :  
at line     122 of function neldermead_variable called by :  
at line       6 of function neldermead_algo called by :  
at line      15 of function neldermead_search called by :  
at line      90 of function fminsearch called by :  
[param,L_Squares,flag,details] = fminsearch(L_Squares,k,options);
at line     136 of exec file called by :    
res_Manip_Couleur-Lumière_Monod_12-12_4.sci', -1

Il est maintenant question d'"illegal input". Je n'y comprend rien...
Merci si vous pouvez m'aider.
Eleasias

[Eleasias]

Mon impression est que l'algorithme de minimisation n'arrive pas à se boucler et pourtant je n'ai pas l'erreur définie dans mon code "le système ne converge pas".
J'ai fait des printf de mon val et en cours de calcul il diminue bien au fil des itérations, puis plante pour on ne sait quelle raison.

Connaissez vous éventuellement un moyen de minimiser en utilisant un autre algorithme sous Scilab?