Mein MATLAB Forum - goMatlab.de

Mein MATLAB Forum

 
Gast > Registrieren       Autologin?   
Bücher:

Fachkräfte:
Entwicklungsingenieur (m/w) für modellbasierte Softwareentwicklung
Modellbasierte Softwareentwicklung mit MATLAB/Simulink und dSPACE TargetLink
Elektronische Fahrwerksysteme GmbH - Ingolstadt

Softwareentwickler (m/w) für automatische Codegenerierung
Softwareentwicklung mit MATLAB/Simulink und dSPACE TargetLink im Bereich Fahrwerkregelsysteme
Elektronische Fahrwerksysteme GmbH - Ingolstadt

Systemingenieur (m/w) Funktionsentwicklung Automotive
Konzeption und Spezifikation von spezifischen Funktionen (Mess- und Regelalgorithmen) in elektronischen Steuergeräten
DRÄXLMAIER Group - Vilsbiburg bei Landshut

Ingenieur (w/m) modellbasierte Entwicklung / Embedded Programmierung
Organisation und Durchführung von Projekten zur Integration von TargetLink in aktuelle Entwicklungsprozesse
dSPACE GmbH - Pfaffenhofen

Embedded Software- Entwickler (m/w) - Hardwarenahe Software für innovative Steuerungssysteme
Design, Implementierung und Test von Embedded Software
DR. JOHANNES HEIDENHAIN GmbH - Traunreut (Raum Rosenheim)

weitere Angebote

Partner:


Vermarktungspartner


Forum
      Option
[Erweitert]
  • Diese Seite per Mail weiterempfehlen
     


Gehe zu:  
Neues Thema eröffnen Neue Antwort erstellen

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

 

Atmer

Gast


Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 24.06.2017, 10:59     Titel: vermutliches zickiger Output-Vektor beim lösen von ODE in SS
  Antworten mit Zitat      
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;
 


Jan S
Moderator

Moderator


Beiträge: 10.481
Anmeldedatum: 08.07.10
Wohnort: Heidelberg
Version: 2009a, 2016b
     Beitrag Verfasst am: 24.06.2017, 15:49     Titel: Re: vermutliches zickiger Output-Vektor beim lösen von ODE i
  Antworten mit Zitat      
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
Private Nachricht senden Benutzer-Profile anzeigen
 
Atmer

Gast


Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 24.06.2017, 19:42     Titel:
  Antworten mit Zitat      
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;

 
 
Jan S
Moderator

Moderator


Beiträge: 10.481
Anmeldedatum: 08.07.10
Wohnort: Heidelberg
Version: 2009a, 2016b
     Beitrag Verfasst am: 26.06.2017, 11:43     Titel:
  Antworten mit Zitat      
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
 

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

Gast


Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 26.06.2017, 13:46     Titel:
  Antworten mit Zitat      
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


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));
 


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.
 
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 goPCB.de


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


Copyright © 2007 - 2017 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.