Mein MATLAB Forum - goMatlab.de

Mein MATLAB Forum

 
Gast > Registrieren       Autologin?   
Bücher:

Fachkräfte:
weitere Angebote

Partner:


Vermarktungspartner


Forum
      Option
[Erweitert]
  • Diese Seite per Mail weiterempfehlen
     


Gehe zu:  
Neues Thema eröffnen Neue Antwort erstellen

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

 

christoph1603
Forum-Newbie

Forum-Newbie


Beiträge: 1
Anmeldedatum: 08.03.19
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 08.03.2019, 15:24     Titel: Zeitäbhängige Variabel in ode-Funktion einfüren
  Antworten mit Zitat      
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;


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
 


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
Private Nachricht senden Benutzer-Profile anzeigen


Neues Thema eröffnen Neue Antwort erstellen



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
.


goMatlab ist ein Teil des goForen-Labels
goForen.de goMATLAB.de goLaTeX.de


 Impressum  | Nutzungsbedingungen  | Datenschutz  | Werbung/Mediadaten | Studentenversion | FAQ | goMatlab RSS Button RSS


Copyright © 2007 - 2019 goMatlab.de | Dies ist keine offizielle Website der Firma The Mathworks
Partner: LabVIEWforum.de

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.