Discussion :

[Eleasias]

Bonjour à tous,
Je cherche à faire coller un modèle biologique à mes données expérimentales à l'aide de Scilab.
Pour être plus précis, j'étudie la croissance de microorganismes en fonction de 2 paramètres (P1 et P2).

(Pour info :
Le type de modèle que je souhaite utiliser est le suivant :
dX/dT = µX * P1/(P1+Ks1)*P2/(P2+Ks2)
Avec X la concentration en biomasse, µ le taux de croissance, P1 et P1 les concentrations de mes paramètres et Ks1 et 2 les constantes de demi saturations utilisées pour limiter ma croissance.)

J'ai 24 expériences avec des couples (P1;P2) différents et mes concentration en microorganismes (8/exp à des temps différents) en cinétique. J'aimerai pouvoir calculer mes 3 constantes, à savoir µ, Ks1 et Ks2. Il semble que la méthode d'optimisation par les moindres carrés soit la plus adaptée pour fiter le modèle.

Je ne sais pas vraiment comment m'y prendre...
En fait pour utiliser la méthode des moindres carrés sur 1 jeu de données, c'est faisable, mais je ne vois pas trop comment faire pour avoir un modèle qui "colle au mieux à l'ensemble des expériences".

Auriez-vous une idée ou de l'expérience en la matière?
Je ne sais pas si je suis assez clair, ni même si vous connaissez les modèles biologiques, alors n'hésitez pas à poser des questions.
En tout cas merci pour votre aide.

Eleasias

[Eleasias]

J'ai tenté d'adapter un code existant à mon problème, mais j'ai quelques problèmes... Pourriez vous m'aider?
Le code existant consistait à calculer 2 coefficients apparaissant dans 2 équations différentielles à l'aide de la méthode des moindres carré et de 2 jeux de données en fonction du temps.

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
// Experimental data
clc
        // data from 25-02-13
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]
 
// Model
function [dy] = mymodel(t,y,flag,u,Ks,DOi,Fluxi)
	dy(1) = u*y(1)*y(2)/(y(2)+Ks)*t 
	dy(2) = =Fluxi*(0.902623-0.02971*y(3)+0.0005421*y(3)*y(3)-0.00000557133*y(3)*y(3)*y(3)+0.0000000299771*y(3)*y(3)*y(3)*y(3) ...
	-0.0000000000653467*y(3)*y(3)*y(3)*y(3)*y(3))
	dy(3) = =12,789+54,7014*DOi+0,000196202*y(1)
	dy = dy'; //transposition
endfunction
// starting guess
u = 0.25;
Ks = 50;
k = [u;Ks];
// --------------------------------------------------------------------------
function [val] = L_Squares(k)
	val=[];
	u = k(1);
	Ks = k(2);
	time = [0 2 4 6 8 10 12 14]';
	y0 = y_exp(1,:);
	t = time;
	y_calc=ode(y0',0,time,mymodel);
	y_calc = y_calc';
	resid = (y_calc - y_exp).*(y_calc - y_exp);
	val = sum(sum(resid));
endfunction
//---------------------------------------------------------------------------
tic
options = optimset("TolFun",1D-02,"TolX",1.D-02,"MaxFunEvals",1.0D+06,"MaxIter",1.0D+06);
[param,L_Squares,flag,details] = fminsearch(L_Squares,k,options);
elapsed_time = toc();
parametros = param'
u = parametros(1); Ks = parametros(2);
clc
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;
time = [0 2 4 6 8 10 12 14]';
t = time;
y=ode([-1 1]',0,time,mymodel);
y = y';
clf
plot(time,y_exp(:,1),"ok")
set(gca(),"auto_clear","off")
plot(time,y_exp(:,2),"ob")
plot(time,y(:,1),"-k")
plot(time,y(:,2),"-Ks")
h1=legend('exp_data_1','exp_data_2','model','model')

Le retour console est le suivant :

  !--error 21 
Index invalide.

at line       4 of function mymodel called by :  
at line       8 of function  called by :  
at line       2 of function  called by :  
at line      59 of function optimbase_function called by :  
at line       6 of function neldermead_costf called by :  
at line       3 of function fun called by :  
at line       6 of function optimsimplex_compsomefv called by :  
at line       5 of function optimsimplex_computefv called by :  
at line      50 of function optimsimplex_pfeffer called by :  
at line      99 of function optimsimplex_new called by :  
at line      31 of function neldermead_startup called by :  
at line       8 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      66 of exec file called by :    
exec('C:\Documents and Settings\...',-1)

Merci si vous acceptez de vous pencher sur mon problème.
Eleasias

[Invité]

Bonjour,

Dans la fonction mymodel, tu as un '=' en trop :

Code :

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

1
2
dy(2) = =Fluxi*(0.902623-0.02971*y(3)+0.0005421*y(3)*y(3)-0.00000557133*y(3)*y(3)*y(3)+0.0000000299771*y(3)*y(3)*y(3)*y(3) ...
	-0.0000000000653467*y(3)*y(3)*y(3)*y(3)*y(3))

Attention aussi que Fluxi est un vecteur de 24 valeurs, le résultat de cette expression sera alors aussi un vecteur d'autant de valeurs, que tu ne pourras pas mettre dans dy(2) qui n'en attend qu'une.