ich möchte die Parameter eines Differenzialgleichungssystems, welches mit einem ODE-solver gelöst wird an Versuchsdaten anpassen. Dazu möchte ich den Levenberg-Marquardt-Algorithmus verwenden. An einem einfacheren Beispiel, ohne ODEs möchte ich mein Problem verdeutlichen.
plot(x,Rexp,'ro'); % plot the experimental data hold on
b0=[11]; % start values for the parameters
options=optimset('Algorithm','levenberg-marquardt');
[b,resnorm,residual,exitflag,output]=lsqnonlin(@recfun,b0,options)% run the lsqnonlin with start value b0, returned parameter values stored in b
Rcal=1./(1+exp(1.0986/b(1)*(x-b(2)))); % calculate the fitted value with parameter b plot(x,Rcal,'b'); % plot the fitted value on the same graph
Rcal=1./(1+exp(1.0986/b(1)*(x-b(2)))); % the calculated value from the model
%y=sum((Rcal-Rexp).^2);
y=Rcal-Rexp;
% the sum of the square of the difference between calculated value and experimental value
laut Doku wird als Syntax vor den Optionen lb und ub erwartet. Wenn du keine unteren und oberen Schranken vorgeben willst, kann man das wohl auf [] setzen. Probier also mal:
ich hatte vorher mit dem Trust-Region-Reflective Algorithmus herumexperimentiert, dem waren die Klammern egal.
Ich habe die Klammern jetzt eingesetzt und den Code incl. Odes und Lsqnonlin gestartet.
Die Fehlermeldung beim Aurfufen erscheint nicht mehr!
MATLAB zeigt jetzt "busy" an (seit 8 Min).
Ich hoffe es wird eine brauchbare Lösung gefunden.
Darf ich dich evtl. später nochmal belästigen, wenn der Solver fertig ist?
schon mal ein Tipp: lsqnonlin versucht mit Differenzenquotienten zu arbeiten. Die Simulationsungenauigkeiten von ODE-Lösern können das ruinieren. Deswegen die Empfehlung,
- mit ODESET die Simulationsgenauigkeit (RELTOL) des ODE-Lösers zu verbessern
und/oder
- mit OPTIMSET die Option DIFFMINCHANGE erhöhen (dann wird dieser Effekt verringert)
Grüße,
Harald
Gast
Gast
Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
Verfasst am: 04.08.2011, 15:20
Titel:
Also was ich gerne mache ist mir die Parameter und die Sum of Squares bei jedem Funktionsaufruf auszugeben.
So kannste zuschauen wie sich
die Werte verändern und das kann je nach Modell recht interessant sein.
danke für den Tipp mit den odeset/optimset Einstellungen. Ich werd das mal ins Auge fassen, wenn der Solver fertig ist. (er ist immer noch nicht fertig )
Als ich die Parameter desselben DGL-Systems mit fminsearch zu optimieren versucht habe, gings erheblich schneller....komisch!
Das mit den Parameter anzeigen lassen muss ich auch mal ausprobieren.
kein Wunder. FMINSEARCH ist "derivative-free", versucht also nicht, Ableitungen anzunähern (und wird vor allem nicht durch unsinnige Ableitungen in die Irre geführt).
Grüße,
Harald
Gast
Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
Verfasst am: 04.08.2011, 17:15
Titel:
Also wenn der jetzt immer noch nicht fertig ist, dann lass dir unbedingt
die Parameter ausgeben um evtl eine Überparametrisierung des
Modells oder so festzustellen. In die Falle rennt so ziemlicher jeder in der Modellierung.
hab jetzt noch ein bischen rumprobiert. Ich glaube, ich übergebe die Daten des ode-Solvers falsch. Ich muss Euch wohl mit meinem Code quälen. Ich hab die unwichtigen Passagen aus lsqparafit1.m rausgelöscht.
Code von lsqparafit1.m
Code:
%%%Parameterfitting mit lsqnonlin%%%
function lsqparafit1
%Festlegung der Zeitschritte für das Lösen des DGL
tspan = HPLC(:,1)'; % Matrix transponieren
%Optionen für die Optimierungsroutine
options = optimset('Algorithm','levenberg-marquardt','DiffMinChange',1e-8,'Display','iter');
%Übergabe der Startparameter an die Least-Square-Routine von lsqnonlin
paropt = lsqnonlin(@minimiz,par0,[],[],options);
%letztmaliges Lösen der DGL mit optimierten Parametern
[t,c_v] = ode45(@lsqparafit2,tspan,c_v0,' ',paropt);
%%% Optimierungsroutine lsqnonlin %%% function lsq=minimiz(par) global c_v0 tspan MAT
%Lösung des DGL-Systems mit Startwerten
[t,c_v] = ode45(@lsqparafit2,tspan,c_v0,' ',par);
%Bewertungskriterium (soll minimiert werden)
tra=c_v';
vekt=[tra(1,:),tra(2,:),tra(3,:)];
lsq = vekt - MAT;
% Ratenkonstanten
kR1= par(2);%kR1=22.3; % (kg/(g*h)) 1/h Reaktionskonstante für r1: Cellulose zu Cellobiose
kR2= par(3);%kR2=7.18; % (kg/(g*h)) 1/h Reaktionskonstante für r2: Cellulose zu Glucose
kR3= par(4);%)285.5; % (1/h) Reaktionskonstante für r3: Cellobiose zu Glucose
% Inhibierung: Glucose (stellv. für C6-Zucker: Glucose,Galaktose, Mannose)
K1Ig= par(5);%0.1; % (g Glucose/kg) Inhibitordissoziationskonstante für r1: Cellulose/Cellobiose
K2Ig= par(6);%0.04; % (g Glucose/kg) Inhibitordissoziationskonstante für r2: Cellulose zu Glucose
K3Ig= par(7);%3.9; % (g Glucose/kg) Inhibitordissoziationskonstante für r3: Cellobiose zu Glucose
% Inhibierung Cellobiose
K1Icb= par(8);%0.015; % (g Cellobiose/kg) Inhibitordissoziationskonstante für r1: Cellulose/Cellobiose
K2Icb= par(9);%132.0; % (g Cellobiose/kg) Inhibitordissoziationskonstante für r2: Cellulose zu Glucose
% Inhibierung: Xylose
Cx=0; % (g Xylose/kg) Xylosekonzentration (stellv. für C5-Zucker: Arabinose,Xylose) in der Lösung
K1Ix=0.1; % (g Xylose/kg) Inhibitordissoziationskonstante für r1: Cellulose/Cellobiose
K2Ix=0.2; % (g Xylose/kg) Inhibitordissoziationskonstante für r2: Cellulose zu Glucose
K3Ix=201.0; % (g Xylose/kg) Inhibitordissoziationskonstante für r3: Cellobiose zu Glucose
% Inhibierung: Lignin (Inhibierung durch Enzymadsorption nach Langmuir)
Cl=20; % (g Lignin/ kg Lösung) Ligninkonzentration 30.4
%lambda=1; % (dimensionslos) Faktor der die bevorzugte Bindung von CBH und EG an Cellulose oder Lignin ausdrückt
KL1L=0.51; % (kg Lösung/g Protein) Langmuir-Verteilungakoeffizient: CHB und EG an Ligninmatrix
KL2L=0.75; % (kg Lösung/g Protein) Langmuir-Verteilungakoeffizient: ß-Glucosidase an Ligninmatrix
Ls1L=0.08607; % (g Protein/g Lignin) maximale Enyzmadsorption CHB und EG je g Lignin
Ls2L=0.1735; % (g Protein/g Lignin) maximale Enyzmadsorption ß-Glucosidase je g Lignin
% Enzymdosierung Cellobiohydrolase CHB und Endoglucanase EG (Hydrolyse zu Cellobiose)
CEs1=45; % (mg Protein/g Cellulose) Enzymdosierung 9 36.59 / 42.1875
Ceg1=CEs1*Cs0/1000; % (g Protein/ kg Lösung) Proteinkonzentration in der Lösung
% Enzymdosierung ß-Glucosidase (Hydrolyse Cellobiose zu Glucose)
CEs2=1; % (mg Protein/g Cellulose) Enzymdosierung 2 8.4 / 2.8125
Ceg2=CEs2*Cs0/1000; % (g Protein/ kg Lösung) Proteinkonzentration in der Lösung
% Substratreaktivität
RS=c_v(1)/Cs0; % (dimensionslos) abnehmende Substratreaktivität mit zunehmendem Abbaugrad der Cellulose
%Adsorptionsteil
% Langmuir Isotherme: CBH und EG an Lignin und Celluloseanteil in der Trockensubstanz
%Langmuirparameter für Adsorption an Cellulose
Ls1=0.04255; % (g Protein/g Cellulose)
KL1=0.6; % (kg Lösung/g Protein)
Ls1t=Ls1+(1-(c_v(1)/Cs0))*(Ls1L-Ls1);
KL1t=KL1+(1-(c_v(1)/Cs0))*(KL1L-KL1);
% Berechnung der freien Konzentration von CBH und EG in der Lösung (g Protein/kg Lösung)
Cefs1=(-(1+Ls1t*KL1t*(c_v(1)+Cl)-Ceg1*KL1t)+((1+Ls1t*KL1t*(c_v(1)+Cl)-Ceg1*KL1t)^2+4*KL1t*Ceg1)^0.5)/(2*KL1t);
Caels=Ceg1-Cefs1;
% Berechnung der freien Konzentration von CBH und EG in der Lösung; Adsorption nur an Ligninanteil (g Protein/kg Lösung)
Cefl1=(-(1+Ls1L*KL1L*Cl-Ceg1*KL1L)+(((1+Ls1L*KL1L*Cl-Ceg1*KL1L)^2)+4*KL1L*Ceg1)^0.5)/(2*KL1L);
Cael3=Ceg1-Cefl1;
% Langmuir Isotherme: CBH und EG an Lignin % Berechnung der gebundenen Konzentration von CBH und EG an Cellulose in der Lösung (g Protein/ kg Lösung)
Ceb1=Caels-Cael3;
% Langmuir Isotherme: ß-Glucosidase an Ligninanteil in der Trockensubstanz % Berechnung der freien Konzentration von ß-Glucosidase in der Lösung (g Protein/kg Lösung)
Cef2=(-(1+Ls2L*KL2L*Cl-Ceg2*KL2L)+((1+Ls2L*KL2L*Cl-Ceg2*KL2L)^2+4*KL2L*Ceg2)^0.5)/(2*KL2L);
% Gleichungen für Reaktionsgeschwindigkeiten, kompetetive Hemmung Glucose, Cellobiose und Xylose (g/(kg*h))
r1=(kR1*Arf*Ceb1*RS*c_v(1))/(1+(c_v(2)/K1Icb)+(c_v(3)/K1Ig)+(Cx/K1Ix));
r2=(kR2*Arf*Ceb1*RS*c_v(1))/(1+(c_v(2)/K2Icb)+(c_v(3)/K2Ig)+(Cx/K2Ix));%Ceb2
r3=(kR3*Arf*Cef2*c_v(2))/(KM*(1+(c_v(3)/K3Ig)+(Cx/K3Ix))+c_v(2));
Der code von lsqparafit2.m ist wahrscheinlich eher unwichtig.
Ich glaube die entscheidende Passage könnte folgende sein, weil ich mir nicht sicher bin was ich an lsqnonlin wie[/color] übergeben muss. Ich übergebe deswegen die Differenz aus 2 Spaltenvektoren. Ein Vektor mit den Messwerten (MAT) und ein Vektor mit den Werten des ode-Solvers (vekt).
Codepassage aus lsqparafit1.m
Code:
%%% Optimierungsroutine lsqnonlin %%% function lsq=minimiz(par) global c_v0 tspan MAT
%Lösung des DGL-Systems mit Startwerten
[t,c_v] = ode45(@lsqparafit2,tspan,c_v0,' ',par);
%Bewertungskriterium (soll minimiert werden)
tra=c_v';
vekt=[tra(1,:),tra(2,:),tra(3,:)];
lsq = vekt - MAT;
die Zeit, die Anwendung von vorne bis hinten durchzugehen, habe ich dann leider doch nicht.
Was du wie übergeben musst, solltest du als Anwender selber am besten wissen. Ich würde aber globale Variablen vermeiden und stattdessen anonyme Function Handles verwenden.
Beispiel: ode45 soll myode aufrufen, die Argumente t, y, param entgegennimmt -->
Der Wert für DiffMinChange ist viel zu niedrig. Dir muss klar sein, dass ode45 nur näherungsweise Lösungen liefert, z.B. mit Toleranz 1e-6 (was ich noch einstellen würde; Default dürfte 1e-3 sein). Ein Differenzenquotient mit Schrittweite 1e-8 ist da unsinnig. Ich würde es mal mit 1e-3 oder 1e-2 für DiffMinChange probieren.
ich werde anonyme Funktionen statt den globalen Parameteren verwenden.
Vielleicht doch noch eine Frage zur Übergabe der Parameter an lsqnonlin, ich werde nicht ganz schlau aus der Anleitung zu lsqnonlin.
Ich übergebe eine Matrix lsq an lsqnonlin.
lsq ist die Differenzmatrix aus vect(der Ergebnismatrix aus dem odesolver) und MAT(den Messwerten). Lsqnonlin soll die Parameter optimieren, so dass diese Differenz minimal wird.
Ist mein Aufruf von lsqnonlin kompatibel mit der Syntax der Funktion,
oder muss ich einen anderen Ausdruck verwenden?
laut Doku sollte die Funktion einen Vektor zurückgeben. Ich weiß nicht, wie sie sich bei Matrizen verhält. Daher würde ich die Matrix von Differenzen in einen Vektor umwandeln, z.B. mit RESHAPE.
Aber warum fragst du?
danke,
das muss ich nächste Woche gleich probieren. Jetzt hab ich leider kein Matlab mehr
Ich kannte die Funtkion reshape bisher nicht.
Ich Frage, weil die Iteration nach dem ersten Schritt nicht weiterläuft. Matlab bricht zwar nicht ab, aber passiert auch nichts Neues (siehe Fehlermeldung oben).
Ich probier den Tips mal aus und geb dann Bescheid.
es liegt aller Wahrscheinlichkeit an DiffMinChange, siehe oben.
Ich würde sicherheitshalber auch die Toleranz von ode45 runtersetzen, siehe auch oben.
Um zu sehen, was MATLAB macht, könntest du dir ja auch in jeder Zielfunktions-Auswertung den verwendeten Parametersatz ins Command Window ausgeben lassen?
also, folgende Tips habe ich umgesetzt (so gut es mir mit meinen beschränkten MATLAB-Kenntnissen möglich ist):
1. mit der Funktion Reshape hab ich Vektoren aus den Matritzen gemacht und dann die Differenz aus den Vektoren an lsqnonlin übergeben
2. ich habe nun mit function handles gearbeitet (sehe ich das richtig, dass ich dann die Variablen immer wieder neu definieren muss?)
3. den Wert für Diffminchange probierte ich als 1e-3 und 1e-2
4. den Wert von Reltol probierte ich als 1e-6 bis 1e-2
Es erscheint leider selbst bei Diffminchange 1e-2 und Reltol 1e-2 eine Meldung:
Du kannst Beiträge in dieses Forum schreiben. Du kannst auf Beiträge in diesem Forum antworten. Du kannst deine Beiträge in diesem Forum nicht bearbeiten. Du kannst deine Beiträge in diesem Forum nicht löschen. Du kannst an Umfragen in diesem Forum nicht mitmachen. Du kannst Dateien in diesem Forum posten Du kannst Dateien in diesem Forum herunterladen
MATLAB, Simulink, Stateflow, Handle Graphics, Real-Time Workshop, SimBiology, SimHydraulics, SimEvents, and xPC TargetBox are registered trademarks and The MathWorks, the L-shaped membrane logo, and Embedded MATLAB are trademarks of The MathWorks, Inc.
)
RSS