Konstantenberechnung über curve fitting

WICHTIG: Der Betrieb von goMatlab.de wird privat finanziert fortgesetzt. - Mehr Infos...

Partner:

Forum

Option

[Erweitert]

• Diese Seite per Mail weiterempfehlen

Gehe zu:

Gast


	Beiträge: ---

	Anmeldedatum: ---

	Wohnort: ---

	Version: ---

Verfasst am: 16.01.2008, 12:49 Titel: Konstantenberechnung über curve fitting

Liebe Matlab Community,
ich habe folgendes Problem und wäre froh, wenn mir einer helfen könnte.

Mein Stoff A reagiert über B zum Produkt C. Die Lösung der Reaktion erfolgt über ein DGL, welches ich mit dem ode45 solver lösen kann. Dafür gebe ich der Funktion die Konstanten k1 und k2 als (Start)werte vor.

Nun möchte ich, dass die Konstanten mithilfe der Least Square Methode (Minimumkriterium: sum( (KonzMesswert - KonzSimulation)^2) ) so angepasst werden, dass die simulierte Lösung möglichst nah an meine Messwerte kommt.

Folgende zwei m-files habe ich geschrieben:
Rkt_ABC (Hauptfile) - hier der Code
clear

%Rkt. A (v13(2,:)) über B (v13(4,:)) zu C (v13(3,:))

%globale Konstanten
global k1 k2

%Datenset
v13=[0 0.0833 0.1667 0.25 0.3333 0.4167 0.5 0.5833 0.6667 0.75 0.9167 1.0833 1.5 2 2.3333 3
6.5 6.19 5.39 4.67 3.88 4.5 3.6 3.08 2.34 1.95 0.93 0.47 0.24 0 0 0
1.88 2.19 2.67 3.35 3.69 4.01 4.51 4.64 5.45 5.44 7.02 8.05 9.5 9.6 9.76 10.02
0.65 0.66 0.81 0.75 0.93 0.9 1.37 1.54 1.78 1.91 1.8 0 0 0 0 0]; %Versuch vom 11.08.2004

%Anfangswerte
t0 = v13(1,1);
c0(1) = v13(2,1);
c0(2) = v13(3,1);
c0(3) = v13(4,1);

%Konstanten per Hand als Startwerte für die spätere Iteration
k1 = 5;
k2 = 1;

%Festlegung der Zeitschritte für das Lösen des DGL
tspan = v13(1,1):.01:v13(1,length(v13));

%Lösung des DGL-Systems mit Startwerten k1 und k2
[t,c] = ode45(@DGL_Rkt_ABC,tspan,c0,[]);

%Anpassen der Konstanten ???
lb = [0,0];
ub = [100,100];
x0 = [k1, k2]; %Startwert für Iteration

% [a]=lsqcurvefit(@function,x0,xdata,ydata,lb,ub);

% k1 = a(1) %in mg(L*h)-1
% k2 = a(2) %in mg(L*h)-1

%plotten
plot(t,c(:,1),'b',v13(1,:),v13(2,:),'bo',t,c(:,2),'r',v13(1,:),v13(4,:),'ro',t,c(:,3),'k',v13(1,:),v13(3,:),'ko');
legend(['A-Simulation'],'A-Experiment',['B-Simulation'],'B-Experiment',['C-Simulation'],'C-Experiment');
xlabel('Zeit in Stunden','Fontsize',10)
ylabel('Konzentrationen in mgL-1','Fontsize',10)
hold on

DGL_Rkt_ABC (Funktion mit dem DGL) - hier der Code
function dcdt = DGL_Rkt_ABC(t,c)

%globale Konstanten
global k1 k2

%DGL-System:

dcdt(1) = - k1*(c(1)); %NH4
dcdt(2) = k1*(c(1))- k2*(c(2)); %NO2
dcdt(3) = k2*(c(2)); %NO3

dcdt=[dcdt(1);dcdt(2);dcdt(3)];

Wie schaffe ich es nun die Konstanten zu optimieren? Brauche ich dazu eine weitere Funktion? Kann ich mit der lsqcurvefit Funktion aus Matlab arbeiten? Wenn ja, wie müssen die Datensätze übergeben werden?

Schon mal im Vorraus Dankend
René

lupie

Forum-Newbie


	Beiträge: 5

	Anmeldedatum: 14.01.08

	Wohnort: Halle(Saale)

	Version: ---

Verfasst am: 20.01.2008, 15:04 Titel:

hi rené,

da du die DGL numerisch berechnest sind einfache fitting-tools in matlab, so glaube ich, nicht geeignet. du hast hier ein optimierungsproblem vorliegen. und dafür gibt es routinen in matlab, wie fminsearch, fmincon...

hab dein programm mal mit einer Optimierung mittels Simplex-Verfahren ausgestattet, welches die fehlerquadratsumme minimiert.
allerdings sehen die verläufe noch nicht wirklich zufriedenstellend aus
dies liegt aber nicht an der optimierung, sondern am mathematischen modell

die Ordnung der Reaktionsgeschw. muss ja nicht immer 1 sein, sie nimmt in der Praxis oft sehr krumme werte an... aber die lassen sich ja auch anpassen Smile

probiers doch mal

hier das programm:

Code:

function main
clear all, close all, clc

global c0 v13 tspan
%Rkt. A (v13(2,:)) über B (v13(4,:)) zu C (v13(3,:))

%Datenset
v13=[0 0.0833 0.1667 0.25 0.3333 0.4167 0.5 0.5833 0.6667 0.75 0.9167 1.0833 1.5 2 2.3333 3
6.5 6.19 5.39 4.67 3.88 4.5 3.6 3.08 2.34 1.95 0.93 0.47 0.24 0 0 0
1.88 2.19 2.67 3.35 3.69 4.01 4.51 4.64 5.45 5.44 7.02 8.05 9.5 9.6 9.76 10.02
0.65 0.66 0.81 0.75 0.93 0.9 1.37 1.54 1.78 1.91 1.8 0 0 0 0 0]; %Versuch vom 11.08.2004

%Anfangswerte
c0(1) = v13(2,1);
c0(2) = v13(4,1);
c0(3) = v13(3,1);

%Konstanten per Hand als Startwerte für die spätere Iteration
k1 = 1.5;
k2 = 4;
par0 = [k1 k2];
%Festlegung der Zeitschritte für das Lösen des DGL
tspan = v13(1,:);
%Optionen für die Optimierungsroutine(optional)
options = optimset('Display','iter','MaxIter',50);
%Übergabe der Startparameter an die Simplex-Optimierungsroutine fminsearch
paropt = fminsearch(@objective,par0,options);
%letztmaliges Lösen der DGL mit optimierten Parametern
[t,c] = ode45(@DGL_Rkt_ABC,tspan,c0,' ',paropt);
%ausgabe der optimierten Parameter
fprintf('k1= %4.2f\nk2= %4.2f\n',paropt(1),paropt(2))
%plotten
plot(t,c(:,1),'b',v13(1,:),v13(2,:),'bo',t,c(:,2),'r',v13(1,:),v13(4,:),'ro',t,c(:,3),'k',v13(1,:),v13(3,:),'ko');
legend(['A-Simulation'],'A-Experiment',['B-Simulation'],'B-Experiment',['C-Simulation'],'C-Experiment');
xlabel('Zeit in Stunden','Fontsize',10)
ylabel('Konzentrationen in mgL-1','Fontsize',10)

%%% Optimierungsroutine fminsearch %%%
function J=objective(par)
global c0 v13 tspan

%Lösung des DGL-Systems mit Startwerten k1 und k2
[t,c] = ode45(@DGL_Rkt_ABC,tspan,c0,' ',par);
%Bewertungskriterium (wird immer minimiert!)aus der Fehlerquadratsumme
J = sum(abs(c(:,1)-v13(2,:)'))+sum(abs(c(:,2)-v13(4,:)'))+sum(abs(c(:,3)-v13(3,:)'));

%%% DGL_Rkt_ABC (Funktion mit dem DGL) - hier der Code %%%
function dcdt = DGL_Rkt_ABC(t,c,par)

k1 = par(1);
k2 = par(2);
%DGL-System:
dcdt(1) = - k1*(c(1)); %NH4
dcdt(2) = k1*(c(1))- k2*(c(2)); %NO2
dcdt(3) = k2*(c(2)); %NO3
dcdt=[dcdt(1);dcdt(2);dcdt(3)];

Funktion ohne Link?

Gast


	Beiträge: ---

	Anmeldedatum: ---

	Wohnort: ---

	Version: ---

Verfasst am: 23.02.2011, 14:26 Titel:

Ist schon ne Weile her, dennoch ein Hinweis zu dem FQS-Term

Code:

J = sum(abs(c(:,1)-v13(2,:)'))+sum(abs(c(:,2)-v13(4,:)'))+sum(abs(c(:,3)-v13(3,:)'));

Funktion ohne Link?

Sollte dieser nicht eigentlich

Code:

J = sum(abs(c(:,1)-v13(2,:)'))+sum(abs(c(:,2)-v13(3,:)'))+sum(abs(c(:,3)-v13(4,:)'));

Funktion ohne Link?

lauten?

Gast


	Beiträge: ---

	Anmeldedatum: ---

	Wohnort: ---

	Version: ---

Verfasst am: 23.02.2011, 14:28 Titel:

Ne ok, hat sich erledigt, habe nicht bemerkt, dass der Datenvektor entsprechend anders ist. Sorry!

Danke @lupie für die schöne Lösung!


Einstellungen und Berechtigungen
Beiträge der letzten Zeit anzeigen:	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

Impressum | Nutzungsbedingungen | Datenschutz | FAQ |

RSS

Hosted by:

Copyright © 2007 - 2024 goMatlab.de | Dies ist keine offizielle Website der Firma The Mathworks

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.