vermutliches zickiger Output-Vektor beim lösen von ODE in S

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

Partner:

Forum

Option

[Erweitert]

• Diese Seite per Mail weiterempfehlen

Gehe zu:

Atmer

Gast


	Beiträge: ---

	Anmeldedatum: ---

	Wohnort: ---

	Version: ---

Verfasst am: 24.06.2017, 11:59 Titel: vermutliches zickiger Output-Vektor beim lösen von ODE in S

Hi,

ich versuche mich gerade an Simulation von Atmung.
Dazu habe mir aus einem Paper Gleichungen gezogen und damit ein state-space-model erstellt, oder eher ich versuche es.

Eines meiner Probleme gerade ist der u-Vektor aus der Gleichung x'=A*x+B*u.
Das will da eben einfach nicht so wie ich es will.
Bisher habe ich mir da dann anderweitig beholfen und habe meinen Input zusammengefasst.
Dadurch habe ich davon natürlich nicht mehr all zu viel, aber ich konnte immerhin etwas plotten.

Noch kurz vorweg etwas zu meinem gleich folgenden Code:
Die Beschriftungen beim Plot sind momentan außer Acht zu lassen,
da ich gerade wild irgendetwas plotte.
Genauer gesagt: Ich plotte momentan das komplette Ergebnis vom ODE-solver, was eben drei Variablen umfasst.
Sollte weiteres Hintergundwissen von Nöten sein, kann ich das gerne geben.

Ich hoffe ihr könnt mir da ein bisschen helfen.

Code:

% Simulation of a multi-compartment model of the lung
% the focus lies on muscle work affecting the change in volume of the lung
clear all
close all

% Parameter definitions
C_AB= 0.1;
C_RC= 0.1;
C_L= 0.2;
C_PL= 0.005;
R_AB= 1.0;
R_RC= 1.0;
R_L= 2.0;

% Input definition: simple rectangle function for muscle contribution
P_AB = @(t)((t>=1)&(t<=2))*0.8;
P_RC = @(t)((t>=1.2)&(t<=2.2))*0.2;
%u = [@(t)((t>=1)&(t<=2))*0.8 ;
% @(t)(((t>=1.2)&(t<=2.2))*0.6)*1];
u=[P_AB;P_RC];

% Definition of the ODE system in standard state-space form
% mathematisch L AB RC
A = [-1/(R_L*C_L)-1/(R_L*C_PL), 1/(R_L*C_PL), 1/(R_L*C_PL);
1/(R_AB*C_PL), -1/(R_AB*C_AB)-1/(R_AB*C_PL), 1/(R_AB*C_PL);
1/(R_RC*C_PL), -1/(R_RC*C_PL), -1/(R_RC*C_RC)-1/(R_RC*C_PL)];

%paper L AB RC
%A = [-1/(R_L*C_L)+1/(R_L*C_PL), -1/(R_L*C_PL), -1/(R_L*C_PL);
% -1/(R_AB*C_PL),-1/(R_AB*C_AB)+1/(R_AB*C_PL),-1/(R_AB*C_PL);
% -1/(R_RC*C_PL),1/(R_RC*C_PL),-1/(R_RC*C_RC)+1/(R_RC*C_PL)];

B1 = [0;
1/(R_AB);
0];

B2 = [0;
0;
1/(R_RC)];

B=[B1,B2];
C = eye(3);
D = 0;
odefun = @(t, Q) (A*Q + B*u(t));

% Standard state space system analysis tools provided by matlab
sys = ss(A, B1, C, D);

% Initial condition
Q0 = [0; 0; 0];

% Time span
tspan = [0; 4];

% Simulate using a numerical ODE solver
[T, Y] = ode45(odefun, tspan, Q0);

% % Plot the results
% figure; %1
% plot(T, [P_AB(T), Y*[0;1;0]]);
% xlabel('Time [s]')
% ylim([-0.1, 1])
% legend('Contribution of abdomen', 'Volume');
% title('Simulation output');
% grid on;
%
%
% figure; %2
% plot(T, [P_RC(T), Y]);
% xlabel('Time [s]')
% ylim([-0.1, 1])
% legend('Contribution of ribcage', 'Volume');
% title('Simulation output');
% grid on;

figure; %3
U =@(t) P_RC(T)+P_AB(T);
plot(T, [U(T),Y]);
xlabel('Time [s]')
ylim([-0.1, 1.3])
legend('sum of muscle contribution','total volume');
title('Simulation output');
grid on;

Funktion ohne Link?

Jan S

Moderator


	Beiträge: 11.057

	Anmeldedatum: 08.07.10

	Wohnort: Heidelberg

	Version: 2009a, 2016b

Verfasst am: 24.06.2017, 16:49 Titel: Re: vermutliches zickiger Output-Vektor beim lösen von ODE

Hallo Atmer,

Was genau ist Deine Frage?

Zitat:

Eines meiner Probleme gerade ist der u-Vektor aus der Gleichung x'=A*x+B*u.
Das will da eben einfach nicht so wie ich es will.

Was willst Du denn?

Gruß, Jan

Atmer

Gast


	Beiträge: ---

	Anmeldedatum: ---

	Wohnort: ---

	Version: ---

Verfasst am: 24.06.2017, 20:42 Titel:

Das Problem mit dem u ist, dass es ein Vektor ist und Matlab, das nicht akzeptier.
Also zumindest nicht auf die Arten, auf die ich es getestet habe.

Meine Gleichung sieht wie folgt aus:
x' = Ax+B1*u1+B2*u2 (A=3x3, x,B1,B2=3x1, u=1x1)
Ich bin mir unsicher, ob das Ergebnis das selbe ist wie bei
x' = Ax+Bu (A=3x3, x,u=3x1, B=3x2)

Außerdem weiß ich nicht genau, wie ich dann jeweils mit der einen oder anderen Variante weiterrechnen soll, vor allem wie ich mit dem B umgehe.

Momentan habe ich eine Version laufen, bei der ich u1 und u2 addiere und B1 und B2 addiere und es dann durch den odesolver jage, das sollte korrekt sein, oder?
Leider bin ich von den Plots nicht so überzeugt.

Code:

% Simulation of a multi-compartment model of the lung
% the focus lies on muscle work affecting the change in volume of the lung

clear all
close all

% Parameter definitions

C_AB= 0.1;
C_RC= 0.1;
C_L= 0.2;
C_PL= 0.005;
R_AB= 1.0;
R_RC= 1.0;
R_L= 2.0;

% Input definition: simple rectangle function for muscle contribution

P_AB = @(t)((t>=1)&(t<=2))*0.8; %u1
P_RC = @(t)((t>=1.2)&(t<=2.2))*0.4; %u2

% Definition of the ODE system in standard state-space form
%x'= Ax + Bu , y=Cx+Du
% system matrix mathematisch L AB RC

A = [-1/(R_L*C_L)-1/(R_L*C_PL), 1/(R_L*C_PL), 1/(R_L*C_PL);
1/(R_AB*C_PL), -1/(R_AB*C_AB)-1/(R_AB*C_PL), 1/(R_AB*C_PL);
1/(R_RC*C_PL), -1/(R_RC*C_PL), -1/(R_RC*C_RC)-1/(R_RC*C_PL)];

%paper L AB RC (erzeugt toggeln)
% A = [-1/(R_L*C_L)+1/(R_L*C_PL), -1/(R_L*C_PL), -1/(R_L*C_PL);
% -1/(R_AB*C_PL),-1/(R_AB*C_AB)+1/(R_AB*C_PL),-1/(R_AB*C_PL);
% -1/(R_RC*C_PL),1/(R_RC*C_PL),-1/(R_RC*C_RC)+1/(R_RC*C_PL)];

%input matrix
B1 = [0;
1/(R_AB);
0];

B2 = [0;
0;
1/(R_RC)];

%B=[B1+B2];
C = eye(3); %wir momentan nicht benutzt
D = 0; %ist immer 0

% Initial condition
Q0 = [0; 0; 0]; %=[x1;x2,x3]=[V_L;V_AB;V_RC]

odefun = @(t, Q) (A*Q + B1*P_AB(t)+B2*P_RC(t));
%odefun = @(t, Q) (A*Q + B*u(t));

% Standard state space system analysis tools provided by matlab
sys = ss(A, B1+B2, C, D);

% Time span
tspan = [0; 4];

% Simulate using a numerical ODE solver
[T, Y] = ode45(odefun, tspan, Q0);

% Plot the results
figure; %1
plot(T, [P_AB(T), Y*[0;1;0]]);
xlabel('Time [s]')
ylim([-0.1, 1 ])
legend('Contribution of abdomen', 'Volume');
title('Simulation output');
grid on;

figure; %2
plot(T, [P_RC(T), Y*[0;0;1]]);
xlabel('Time [s]')
ylim([-0.1, 1])
legend('Contribution of ribcage', 'Volume');
title('Simulation output');
grid on;

figure; %3
plot(T, [P_AB(T)+P_RC(T),Y*[1;0;0]]);
xlabel('Time [s]')
ylim([-0.1, 1.5])
legend('sum of muscle contribution','total volume');
title('Simulation output');
grid on;

Funktion ohne Link?

Jan S

Moderator


	Beiträge: 11.057

	Anmeldedatum: 08.07.10

	Wohnort: Heidelberg

	Version: 2009a, 2016b

Verfasst am: 26.06.2017, 12:43 Titel:

Hallo Atmer,

Denke daran, dass die Leser im Forum nicht die geringste Ahnung davon haben, was Du tust und was damit beabsichtigt ist.

Zitat:

Das Problem mit dem u ist, dass es ein Vektor ist und Matlab, das nicht akzeptiert.

Dazu kann man kaum etwas antworten. Was heißt denn "nicht akzeptiert"? Bekommst Du eine Fehlermeldung? Dann poste sie und füge den relevanten code ein.

Zitat:

Also zumindest nicht auf die Arten, auf die ich es getestet habe.

Welche Arten sind das?

Zitat:

x' = Ax+B1*u1+B2*u2 (A=3x3, x,B1,B2=3x1, u=1x1)
Ich bin mir unsicher, ob das Ergebnis das selbe ist wie bei
x' = Ax+Bu (A=3x3, x,u=3x1, B=3x2)

Das kommt auch darauf an, was "u1" und "u2" bedeutet, wenn "u=1x1".

Zitat:

Außerdem weiß ich nicht genau, wie ich dann jeweils mit der einen oder anderen Variante weiterrechnen soll, vor allem wie ich mit dem B umgehe.

Das lässt sich auch nicht beantworten. Konkrete Fragen sind besser.

Zitat:

Momentan habe ich eine Version laufen, bei der ich u1 und u2 addiere und B1 und B2 addiere und es dann durch den odesolver jage, das sollte korrekt sein, oder?

Ob es korrekt ist, kannst nur Du entscheiden. Solange nicht klar ist, was "u1" und "u2" ist, kann man das als Leser kaum beurteilen.

Zitat:

Code:

P_AB = @(t)((t>=1)&(t<=2))*0.8; %u1
P_RC = @(t)((t>=1.2)&(t<=2.2))*0.4; %u2

Funktion ohne Link?

Das sieht nach unstetigen Funktionen aus. ODE45 kann nur stetig differnzierbare Funktionen zuverlässig integrieren. Jede Sprung sorgt für ein chotishes Verhalten der Schrittweiten-Kontrolle: Entweder stoppt die Integration mit einer Fehlermeldung, wenn Du Glück hast. Oder Du bekommst ein unzuverlässiges Ergebnis ohne jeden Hinweis.
Als wissenschaftliches Endergebnis kann man dann eine solche Integration nicht mehr betrachten. Als Betreuer einer Doktorarbeit würde ich das einem Studenten "um die Ohren hauen".

Ich finde anonyme Funktionen immer deutlich unübersichtlicher für umfangreiche Formeln und bevorzuge herkömmliche Funktionen. Dann ließen sich die Unstetigkeiten einfach mit einer Schleife über die Teil-Intervalle behandeln. Du könntest den Debugger verwenden um zeilenweise durch den Code zu gehen. Wenn Du dann die Fehlermeldungen zu den entsprechenden Zeilen postest, ist eine konkrete Hilfe möglich.

Gruß, Jan

Atmer

Gast


	Beiträge: ---

	Anmeldedatum: ---

	Wohnort: ---

	Version: ---

Verfasst am: 26.06.2017, 14:46 Titel:

Hallo Jan,

ich saß einfach zu lange vor diesem Code, dass ich eigentlich selbst gar nicht mehr genau wusste, was ich tun wollte und welcher Teil nun überhaupt auf welchem Stand ist, weil ich überall wild umher geändert und nicht genügend Zwischenversionen gespeichert habe.

Zitat:

Denke daran, dass die Leser im Forum nicht die geringste Ahnung davon haben, was Du tust und was damit beabsichtigt ist.

Mein Versuch war die Atmung zu simulieren, mit Matlab, mit einem System von Differentialgleichungen, mit Hilfe eines odesolvers.
Dafür benötigte ich zunächst einen Zustandsraum der Form x'=Ax+Bu.
A - Systemmatrix
x - Zustände
B - Eingangsmatrix
u - Eingänge

Zitat:

Das kommt auch darauf an, was "u1" und "u2" bedeutet, wenn "u=1x1".

Hier habe ich mich verschrieben, u1 und u2 sind jeweils 1x1

Zitat:

Momentan habe ich eine Version laufen, bei der ich u1 und u2 addiere und B1 und B2 addiere und es dann durch den odesolver jage, das sollte korrekt sein, oder?

An dieser Stelle wollte ich wissen, ob es legitim ist die B-Matrix aufzuteilen,
denn rein mathematisch ist es vollkommen unerheblich, Hauptsache die Dimension stimmt eben im Nachhinein immernoch.

Code:

P_AB = @(t)((t>=1)&(t<=2))*0.8; %u1
P_RC = @(t)((t>=1.2)&(t<=2.2))*0.4; %u2

Funktion ohne Link?

Das sind Rechteckfunktionen und die sind meines Wissens nach nicht stetig.
Tatsächlich funktioniert das eigentlich trotzdem sehr gut.
Ich habe zumindest ein Beispiel-Code, in dem es um Medikamentenaufnahme in Körperregionen geht und da wird auch nur kurz ein Rechtecksignal reingegeben, um zu bedeuten, dass das Medikament in den Körper gelangt ist (bspw durch orale Einnahme).

Mit meinem Betreuer zusammen habe ich mir für die Simulation aber genau diesen Eingang zusammen überlegt.
Vielleicht sollte ich aber fragen, ob wir daraus lieber eine Gauß-kurve oder Ähnliches machen. (Ich hoffe ich meine da gerade die Gaußglocke)

Ein großes B habe ich übrigens mittlerweile hinbekommen,
bzw habe mir da helfen lassen.

Code:

P = @(t) [P_AB(t); P_RC(t)];
B = [B1, B2];
odefun = @(t, Q) (A*Q + B*P(t));

Funktion ohne Link?

Zur Erinnerung:
P_AB und P_RC waren beides Funktionen abhängig von der Zeit und sollten meine Eingangsgrößen sein.
Ich konnte sie allerdings nicht gemeinsam in einen Eingangsvektor bringen.
Das musste über eine weitere Funktion geschehen, wie man im Codefragment hier sieht.


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.