Optimierung der Anfangswerte einer DGL

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

Partner:

Forum

Option

[Erweitert]

• Diese Seite per Mail weiterempfehlen

Gehe zu:

ElArminio

Forum-Newbie


	Beiträge: 9

	Anmeldedatum: 19.05.10

	Wohnort: ---

	Version: ---

Verfasst am: 14.06.2010, 00:20 Titel: Optimierung der Anfangswerte einer DGL

Hallo!

Ich habe folgendes Problem:
Ich möchte mit der Funktion fminsearchbnd die Anfangswerte meines DGL-Systems (chem. Folgereaktion) so optimieren, dass die Selektivität 1 (siehe Code) wird.

Code:

function Rahmen_selekt

global tspan v Vr

c0_Start=[0.456 0 0];

%Systembedingungen

k_opt=[0.04977 0.049833];
Vr=10; % Reaktorvolumen [m³]
v=0.1; % Volumenstrom/Durchsatz [m³/s]

LB=[0 0 0]; %Unter- und...
UB=[2 0 0]; %...Obergrenze der Anfangskonzentration c_A und B und C

tspan=0:0.01:200;

%Optionen für die Optimierungsroutine
options = optimset('Display','iter','MaxIter',1000);

c0_neu = fminsearchbnd(@Selekt,c0_Start,LB,UB,options);

c_A_opt=c0_neu(1);

fprintf('Optimale Eintrittskonzentration von A: \ncA_opt= %4.3f\n\n', c_A_opt)

[T,C]=ode45(@system,tspan,c0_neu,k_opt,v,Vr);

plot(T,C)

function S_Cinv = Selekt(par)
global tspan k_opt

[T,C]=ode45(@system,tspan,par,' ',k_opt);
S_C=C(tspan(end),3)/(c0(1)-C(tspan(end),1)); %Selektivität liegt zwischen 0 und 1, soll möglichst 1 werden
S_Cinv = 1-S_C; %"Umwandlung", damit minimiert wird

function dc = system(t,c,par,v,Vr)
global k_opt c0

v=0.1; Vr=10;
k_opt=[0.04977 0.049833];

c0(1) = par(1);
c0(2) = par(2);
c0(3) = par(3);

%DGL-System:
dc=[-k_opt(1)*c(1)+v/Vr*(c0(1)-c(1));
-k_opt(2)*c(2)+k_opt(1)*c(1)+v/Vr*(c0(2)-c(2));
k_opt(2)*c(2)+ v/Vr*(c0(3)-c(3))];

Funktion ohne Link?

Als Fehlermeldung bekomme ich:

Code:

??? Attempted to access par(3); index out of bounds because numel(par)=2.

Funktion ohne Link?

Ich verstehe nicht, an welcher Stelle die Anzahl der Komponenten (numel) festgelegt wird. Sieht da jemand durch und kann weiterhelfen?

Oder anders gefragt: Wie würde man in so einem Fall vorgehen mit der Fehlersuche?

Vielen Dank!
Armin
_________________

Und das letzte was wir hören werden, wenn die Welt explodiert sind die Worte eines Experten, der sagt: "Das ist technisch unmöglich."

Thomas84

Forum-Meister


	Beiträge: 546

	Anmeldedatum: 10.02.10

	Wohnort: ---

	Version: ---

Verfasst am: 14.06.2010, 06:56 Titel:

ich glaub der ode45 Aufruf ist nicht ganz richtig. Ich schreib immer sowas:

Code:

[T,C]=ode45(@(t,c) system(t,c,c0_neu,v,Vr),tspan);

Funktion ohne Link?

ElArminio

Themenstarter

Forum-Newbie


	Beiträge: 9

	Anmeldedatum: 19.05.10

	Wohnort: ---

	Version: ---

Verfasst am: 14.06.2010, 07:28 Titel:

Also das ist mir neu...
In der Befehlsdoku von Mathworks steht meine Variante auch drin. Hat auch bisher so hingehauen.

Ich habe auch bereits eine ähnliche Optimierungsaufgabe fertig, bei der es darum geht, den zusätzlichen Parameter k mit Messwerten anzupassen (ergo Minimierung der Fehlerquadratsumme).
[/url]
_________________

Und das letzte was wir hören werden, wenn die Welt explodiert sind die Worte eines Experten, der sagt: "Das ist technisch unmöglich."

Thomas84

Forum-Meister


	Beiträge: 546

	Anmeldedatum: 10.02.10

	Wohnort: ---

	Version: ---

Verfasst am: 14.06.2010, 07:38 Titel:

an welcher Stelle denn? Im Beispiel 3 stehts z.B.:

Code:

[T Y] = ode45(@(t,y) myode(t,y,ft,f,gt,g),Tspan,IC); % Solve ODE

Funktion ohne Link?

ElArminio

Themenstarter

Forum-Newbie


	Beiträge: 9

	Anmeldedatum: 19.05.10

	Wohnort: ---

	Version: ---

Verfasst am: 14.06.2010, 18:37 Titel:

Dein Syntax geht auch, ein Beispiel für "meine" Art des Aurufs steht im ersten Beispiel.

Code:

[T,Y] = ode45(@rigid,[0 12],[0 1 1],options);

Funktion ohne Link?

Allerdings hat Matlab denke ich schon vorher ein Problem, bevor es den ode45-Solver bemüht.
Hat jemand eine Idee?

Gruß Armin
_________________

Und das letzte was wir hören werden, wenn die Welt explodiert sind die Worte eines Experten, der sagt: "Das ist technisch unmöglich."

ElArminio

Themenstarter

Forum-Newbie


	Beiträge: 9

	Anmeldedatum: 19.05.10

	Wohnort: ---

	Version: ---

Verfasst am: 14.06.2010, 18:40 Titel:

Dein Syntax geht auch, ein Beispiel für "meine" Art des Aurufs steht im ersten Beispiel.

Code:

[T,Y] = ode45(@rigid,[0 12],[0 1 1],options);

Funktion ohne Link?

Thomas84

Forum-Meister


	Beiträge: 546

	Anmeldedatum: 10.02.10

	Wohnort: ---

	Version: ---

Verfasst am: 15.06.2010, 06:32 Titel:

Das ist genau der gleiche Syntax:

y = ode45(funktion,zeitvektor,anfangswerte,optionen)

Der Fehler tritt zumindest auf weil par nicht richtig übergeben wird. Versteh sowieso nicht warum du v und Vr mit übergibst und dann überschreibst. Außerdem kann man auch die ganzen globalen Variablendefinitionen weglassen. Dann wird es etwas übersichtlicher.

ElArminio

Themenstarter

Forum-Newbie


	Beiträge: 9

	Anmeldedatum: 19.05.10

	Wohnort: ---

	Version: ---

Verfasst am: 15.06.2010, 23:22 Titel:

Heureka, es funktioniert! Very Happy

Es lag wirklich an der falschen Übergabe von und dann noch ein paar Kleinigkeiten. Ein paar globals brauche ich aber trotzdem zB für tspan.
Vielen Dank thomas84 für deinen Rat!

Für alle, die es interessiert, ist hier der Code nochmal:

Code:

function Rahmen_selekt

%Systembedingungen

c0_Start=[0.456 0 0];
Vr=10; % Reaktorvolumen [m³]
v=0.1; % Volumenstrom/Durchsatz [m³/s]

LB=[0 0 0]; %Unter- und...
UB=[2 0 0]; %...Obergrenze der Anfangskonzentration c_A und B und C

tspan=0:0.01:200;

%Optionen für die Optimierungsroutine
options = optimset('Display','iter','MaxIter',1000);

%Übergabe der Startparameter an die Simplex-Optimierungsroutine fminsearchbnd

c0_neu = fminsearchbnd(@Selekt,c0_Start,LB,UB,options);

c_A_opt=c0_neu(1);

fprintf('Optimale Eintrittskonzentration von A: \ncA_opt= %4.3f\n\n', c_A_opt)

[T,C]=ode45(@(t,c) system(t,c,c0_neu,v,Vr),tspan,c0_neu);

S_Copt=C(tspan(end),3)/(c0_Start(1)-C(tspan(end),1));
fprintf('Optimale Selektivität von C: \nS_Copt= %4.3f\n\n', S_Copt)

figure
title('Lösung des DGL-Systems mit optimierter Eintrittskonzentration')
plot(T,C)
legend('c_A','c_B ','c_C');
xlabel('Zeit t [s]');
ylabel('Konzentration [kmol/m³]');
text(3e1,0.47,['$c_{A0,opt}=$',num2str(c_A_opt),' $\frac{kmol}{m^3}$'],'interpreter','latex','FontSize',13);

function S_Cinv = Selekt(par)
global tspan v Vr

c0_Start=[0.456 0 0];
[T,C]=ode45(@(t,c) system(t,c,par,v,Vr),tspan,par);
S_C=C(tspan(end),3)/(c0_Start(1)-C(tspan(end),1));
S_Cinv = (1-S_C)^2;

function dc = system(t,c,par,v,Vr)

k_opt=[0.04977 0.049833]; % [1/s] manuelle Übergabe aus Rahmenskript

c0(1) = par(1);
c0(2) = par(2);
c0(3) = par(3);

%DGL-System:
dc=[-k_opt(1)*c(1)+v/Vr*(c0(1)-c(1));
-k_opt(2)*c(2)+k_opt(1)*c(1)+v/Vr*(c0(2)-c(2));
k_opt(2)*c(2)+ v/Vr*(c0(3)-c(3))];

Funktion ohne Link?

Gruß Armin
_________________

Und das letzte was wir hören werden, wenn die Welt explodiert sind die Worte eines Experten, der sagt: "Das ist technisch unmöglich."


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 - 2025 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.