Zeitäbhängige Variabel in ode-Funktion einfüren

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

Partner:

Forum

Option

[Erweitert]

• Diese Seite per Mail weiterempfehlen

Gehe zu:

christoph1603

Forum-Newbie


	Beiträge: 1

	Anmeldedatum: 08.03.19

	Wohnort: ---

	Version: ---

Verfasst am: 08.03.2019, 15:24 Titel: Zeitäbhängige Variabel in ode-Funktion einfüren

Hallo Zusammen,
ich versuche gerade seit einer Woche ein Problem in diesem Programm zu lösen.

Bevor ich zu meinen Problem komme, etwas allgemeines zum Programm.
Das Programm diehnt zur Schätzung von Parametern einer Gleichung. Dazu werden Daten (Konz und Zeit) gemessen, implementiert und eine nicht lineare regression durchgeführt, an der die unbekannten Parameter einer Gleichung angepasst werden.

Ich möchte jetzt die verwendeten Daten noch um die Daten in VR erweitern. VR soll dabei für jeden Zeitstempel in der Gleichung eingesetzt werden. VR und Zeit haben als vektor identische größe.

Leider funktioniert meine Idee die ich hatte nicht. Kann mir wer einen Tipp geben wie ich das meistern kann? Ich wäre euch sehr verbunden.

Code:

clear
clear global;
close;
clear;
global Konz TempC pbar yest
global nBeg nEnd n_V VR
% Anzahl der Versuche
n_V=16;
% Anfangskonzentration [mol/m³]
ca0 =1429.526476;
% Probenameraster [min]
t_V1=[0;5;10;15;30;45;60;75;105;135];
t_V2=[0;5;10;15;30;45;60;75;85;87.5000000000000];
t_V3=[0;5;10;15;30;45;50;55];
t_V4=[0;5;10;15;30;35;40];
t_V5=[0;5;10;15;25;40;50;52;55;57];
t_V6=[0;5;10;15;25;35;40;42];
t_V7=[0;5;10;15;20;25;30;32];
t_V8=[0;5;10;15;20;22;25];
t_V9=[0;5;10;15;20;25;30;32.5000000000000;35];
t_V10=[0;5;10;15;20;22.5000000000000;25];
t_V11=[0;5;10;15;17.5000000000000;20;22.5000000000000];
t_V12=[0;5;10;12.5000000000000;15];
t_V13=[0;5;10;15;20;22.5000000000000;25;27.5000000000000];
t_V14=[0;5;10;15;17.5000000000000;20];
t_V15=[0;5;10;12.5000000000000;15];
t_V16=[0;2.50000000000000;5;7.50000000000000;10;12.5000000000000];
% Versuchstemperaturen [°C]
TempC=[30 40 50 60 30 40 50 60 30 40 50 60 30 40 50 60];
%Versuchsdruck in bar
pbar=[2 2 2 2 4 4 4 4 6 6 6 6 8 8 8 8];
%Gemessene Konzentrationen in den Proben [mol/m³]
Konz_V1=[1429.52647600000;1272.86936600000;1132.14553700000;1022.86846600000;851.455537200000;695.951549900000;556.262584000000;432.304509400000;231.292716600000;30.7992670400000];
Konz_V2=[1429.52647600000;1210.29419100000;1100.90792800000;976.082167600000;727.027560200000;540.757032000000;354.930811000000;154.118618600000;30.7992670400000;0];
Konz_V3=[1429.52647600000;1210.29419100000;1007.26992100000;867.022992400000;556.262584000000;246.736733900000;107.850917700000;15.3981063400000];
Konz_V4=[1429.52647600000;1132.14553700000;913.743991000000;742.570209700000;308.543492500000;138.692990300000;15.3981063400000];
Konz_V5=[1429.52647600000;1163.39563000000;991.674487900000;835.891186000000;618.315665000000;339.465298300000;154.118618600000;77.0210843700000;46.2034830000000;0];
Konz_V6=[1429.52647600000;1163.39563000000;976.082167600000;773.664802000000;509.755183500000;231.292716600000;77.0210843700000;15.3981063400000];
Konz_V7=[1429.52647600000;1132.14553700000;867.022992400000;649.360742300000;447.788482500000;262.183819300000;61.6107551400000;15.3981063400000];
Konz_V8=[1429.52647600000;1085.29380100000;758.115956700000;494.258885100000;231.292716600000;92.4344715900000;15.3981063400000];
Konz_V9=[1429.52647600000;1100.90792800000;867.022992400000;649.360742300000;447.788482500000;308.543492500000;138.692990300000;77.0210843700000;15.3981063400000];
Konz_V10=[1429.52647600000;1054.07490200000;727.027560200000;447.788482500000;231.292716600000;107.850917700000;15.3981063400000];
Konz_V11=[1429.52647600000;1085.29380100000;649.360742300000;354.930811000000;184.979063800000;46.2034830000000;0];
Konz_V12=[1429.52647600000;898.167218400000;354.930811000000;123.270423600000;0];
Konz_V13=[1429.52647600000;991.674487900000;727.027560200000;432.304509400000;215.851766500000;107.850917700000;46.2034830000000;0];
Konz_V14=[1429.52647600000;929.323871300000;587.282948100000;215.851766500000;77.0210843700000;0];
Konz_V15=[1429.52647600000;898.167218400000;354.930811000000;107.850917700000;0];
Konz_V16=[1429.52647600000;1038.47012700000;649.360742300000;370.399398200000;123.270423600000;0];
% Reaktorvolumen während der Versuche [m³]
V_V1=[0.000200000000000000;0.000199999999999999;0.000194300000000000;0.000187600000000000;0.000180900000000000;0.000174100000000000;0.000167500000000000;0.000160800000000000;0.000154200000000000;0.000147300000000000];
V_V2=[0.000200000000000000;0.000199999999999999;0.000193500000000000;0.000186800000000000;0.000179900000000000;0.000173200000000000;0.000166300000000000;0.000159700000000000;0.000152700000000000;0.000145900000000000];
V_V3=[0.000200000000000000;0.000199999999999999;0.000193200000000000;0.000186300000000000;0.000179200000000000;0.000172500000000000;0.000165800000000000;0.000158900000000000];
V_V4=[0.000200000000000000;0.000199999999999999;0.000192900000000000;0.000186000000000000;0.000179000000000000;0.000172200000000000;0.000165200000000000];
V_V5=[0.000200000000000000;0.000199999999999999;0.000193500000000000;0.000186800000000000;0.000180100000000000;0.000173400000000000;0.000166800000000000;0.000160300000000000;0.000153400000000000;0.000146300000000000];
V_V6=[0.000200000000000000;0.000199999999999999;0.000193200000000000;0.000186600000000000;0.000180000000000000;0.000173200000000000;0.000166200000000000;0.000159100000000000];
V_V7=[0.000200000000000000;0.000199999999999999;0.000193100000000000;0.000186400000000000;0.000179800000000000;0.000172900000000000;0.000166400000000000;0.000159400000000000];
V_V8=[0.000200000000000000;0.000199999999999999;0.000193200000000000;0.000186400000000000;0.000179600000000000;0.000173000000000000;0.000164700000000000];
V_V9=[0.000200000000000000;0.000199999999999999;0.000193300000000000;0.000186700000000000;0.000179900000000000;0.000173200000000000;0.000166500000000000;0.000159600000000000;0.000152900000000000];
V_V10=[0.000200000000000000;0.000199999999999999;0.000193300000000000;0.000186700000000000;0.000180000000000000;0.000173400000000000;0.000166500000000000];
V_V11=[0.000200000000000000;0.000199999999999999;0.000193500000000000;0.000186600000000000;0.000179800000000000;0.000173100000000000;0.000166200000000000];
V_V12=[0.000200000000000000;0.000199999999999999;0.000193200000000000;0.000186600000000000;0.000179800000000000];
V_V13=[0.000200000000000000;0.000199999999999999;0.000193500000000000;0.000186900000000000;0.000180400000000000;0.000173900000000000;0.000167400000000000;0.000160900000000000];
V_V14=[0.000200000000000000;0.000199999999999999;0.000193300000000000;0.000186900000000000;0.000180400000000000;0.000173700000000000];
V_V15=[0.000200000000000000;0.000199999999999999;0.000193300000000000;0.000186600000000000;0.000180000000000000];
V_V16=[0.000200000000000000;0.000199999999999999;0.000193400000000000;0.000186600000000000;0.000179800000000000;0.000173000000000000];
%Zusammenfassen der Daten
Zeit=vertcat(t_V1, t_V2, t_V3, t_V4, t_V5, t_V6, t_V7, t_V8, t_V9, t_V10, t_V11, t_V12, t_V13, t_V14, t_V15, t_V16);
Konz=vertcat(Konz_V1, Konz_V2, Konz_V3, Konz_V4, Konz_V5, Konz_V6, Konz_V7, Konz_V8, Konz_V9, Konz_V10, Konz_V11, Konz_V12, Konz_V13, Konz_V14, Konz_V15, Konz_V16);
VR=vertcat(V_V1, V_V2, V_V3, V_V4, V_V5, V_V6, V_V7, V_V8, V_V9, V_V10, V_V11, V_V12, V_V13, V_V14, V_V15, V_V16);
% Initialisierung Iteration
nBeg=[1 11 21 29 36 46 54 62 69 78 85 92 97 105 111 116];
nEnd=[10 20 28 35 45 53 61 68 77 84 91 96 104 110 115 121];
%Startwerte
%Ansatz 1
k1= 20000;
Ea1= 30000;
kAu= 0.0001;
kHu= 0.0001;
dA=-10000;
dH=-10000;
% Ober/Untere Grenzen
%Ansatz 1
ub=[50000 60000 0.01 0.01 -2000 -2000];
lb=[10000 20000 0.000001 0.000001 -30000 -30000];
% Option für lsqcurvefit
options=optimset;
options = optimset(options,'Display','iter-detailed','TolFun',1e-15,'TolX',4e-6); %Ansatz 1
options.MaxFunEvals = 2000;
% Parameterschätzung
[par,resnorm,residual,exitflag,output,lambda,jacobian] = lsqcurvefit(@ParSchn, [k1 Ea1 kAu kHu dA dH], Zeit, Konz, VR, lb,ub,options);
%Konfidenzintervall der Parameter
ci=nlparci(par,residual,'jacobian',jacobian);
%Ausgabe Umsatz
Uest=(ca0-yest(:,1))/ca0;

Funktion ohne Link?

Code:

function [cA]=STRn(t,c,VR)
global k1 Ea1 Temp kAu kHu dA dH p
mcat=0.5; % Katalysatormasse [g]
%Löslichkeit H2 [mol/m³]
cH2=p*(0.0145*Temp-1.6985);
% Arrhenius Gleichung
k=k1*exp(-Ea1/(8.3143*Temp));
kA=kAu*exp(-dA/(8.3143*Temp));
kH=kHu*exp(-dH/(8.3143*Temp));
% Reaktionsgeschwindigkeit
%Ansatz 1
r=(mcat/VR(1) * k * abs(c(1))*kA*kH*cH2)/((1+abs(c(1))*kA)*(1+kH*cH2));
% Stoffbilanz STR
cA=-r;
end

Funktion ohne Link?

Code:

function yMod=ParSchn(par,Zeit)
global k1 Ea1 Temp kAu kHu dA dH p yest
global Konz TempC pbar VR
global nBeg nEnd n_V
k1= par(1);
Ea1= par(2);
kAu= par(3);
kHu= par(4);
dA=par(5);
dH=par(6);
% Wdh. für alle Versuche
for i=1:n_V
Temp=TempC(i)+273.15; % Berechnung T in K
p=pbar(i);
% Aufruf Lösung DGL für STR
[x,y]=ode23s(@STRn,Zeit(nBeg(i):nEnd(i)),Konz(nBeg(i)), VR(nBeg(i):nEnd(i)));
figure(i);plot(x,y(:,1)); % Plot Regression
hold on; % Figure für weitere Daten offen halten
plot(Zeit(nBeg(i):nEnd(i)),Konz(nBeg(i):nEnd(i),1),'o'); %Plot Messdaten
% Titel und Achsenbeschriftung einfügen
title(['Temperatur: ' num2str(Temp) 'K ','Druck: ' num2str(p) 'bar ']);
xlabel('t [min]')
ylabel('c (mol/m³)')
hold off;
% Model Vorhersagen werden in Matrix yest gespeichert
yMod(nBeg(i):nEnd(i),1:1)=y(:,1);
yest=yMod;
end
end

Funktion ohne Link?


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.