Anpassung von Parametern von ODEs bei Verwendung von ODE-Löser octave/matlab

Question

Ich verwende OdePkg in Octave, um ein System von steifen ODEs, z. durch ode5r:

function yprime = myODEs(t,Y,param) 
    yprime = [ 
     - param(1) * Y(1);      # ODE for Y(1) 
     param(1) * Y(1) - param(2) Y(2) * Y(3); # ODE for Y(2) 
     param(2) Y(2) * Y(3)     # ODE for Y(3) 
               # etc. 
]; 

time_span = [1, 24]   # time span of interest 
Y0  = [1.0, 1.1, 1.3] # initial values for dependent variables Y 
param  = [7.2, 8.6, 9.5] # parameters, to be optimized 

[t, Y] = ode5r(@myODEs, time_span, Y0, ..., param); 

Der Solver speichert die abhängigen Variablen Y in einer Matrix mit Bezug auf die Zeit t (Vektor):

t  Y(1) Y(2) Y(3) 
0.0 1.0 1.1 1.3 
0.1 ... ... ... 
0.5 ... ... ... 
0.9 ... ... ... 
... ... ... ... 
4.0 ... ... ... 
... ... ... ... 
24.0 ... ... ... 

I die Parameter in param passen soll, so dass die resultierenden Variablen Y besten meine Referenzwerte passen, zB:

t   Y(1) Y(2) Y(3) 
0.5  1.1 N/A N/A 
1.0  1.9 N/A N/A 
4.0  2.3 2.7 2.1 
5.0  N/A 2.6 2.2 
24.0  0.9 1.5 2.0 

Welche Octave/Matlab (andere Sprachen sind willkommen) Routine kann führen Sie einen Multi-Parameter (Least-square/Spline) passen? Wie ist es möglich Parametersätze für verschiedene Anfangswerte zu kombinieren Y0 in der Anpassung? Ich würde mich freuen, wenn Sie mir ein paar Tipps und Möglichkeiten geben könnten.

Mit freundlichen Grüßen, Simon

Answer 1

Dies sollte relativ einfach mit scipy sein. scipy.optimize.leastsq() nimmt eine Funktion, die ein Array von Residuen für einen gegebenen Parametervektor zurückgeben soll. Es minimiert die Summe der Quadrate der Residuen. Um mehrere Datensätze mit unterschiedlichen Anfangswerten zu verarbeiten, führen Sie die ODE nur einmal für jedes Dataset aus, berechnen Sie die Residuen für jedes Dataset/Ausführungspaar und verketten Sie dann die restlichen Vektoren miteinander. Hier ist eine grobe Skizze:

import numpy 
from scipy import integrate, optimize 

# The initial guess. 
p0 = numpy.array([7.2, 8.6, 9.5]) 

# The collected datasets. 
# A list of (t, y0, y) tuples. 
# The y's are all (len(y0), len(t))-shaped arrays. The output of 
# integrate.odeint is also in this format. 
datasets = [...] 

def odes(y, t, params): 
    dydt = [ 
     -params[0] * y[0], 
     params[0]*y[0] - params[1]*y[1]*y[2], 
     params[1]*y[1]*y[2], 
    ] 
    return np.array(dydt) 

def residuals(params, datasets): 
    res = [] 
    for t, y0, y in datasets: 
     res.append(integrate.odeint(odes, y0, t, args=(params,)) - y) 

    # Stack them horizontally and flatten the array into the expected vector. 
    # You're on your own for handling missing data. Look into the numpy.ma 
    # module. 
    all_residuals = numpy.hstack(res).ravel() 
    return all_residuals 

opt_params, err = optimize.leastsq(residuals, p0, args=(datasets,)) 

Answer 2

Meinen Sie, dass jede Funktion y (t) ausgestattet werden muss? In diesem Fall funktioniert ein Leasingquadrat oder eine Spline-Anpassung für jeden Yi-Satz gegen die Zeit gut. Kann nicht sagen, welches wäre am besten, ohne Ihre Daten zu sehen.

Sie müssen sich eine andere unabhängige Variable einfallen lassen, wenn Sie eine Kurve über alle Werte von Yi für einen bestimmten Zeitpunkt anpassen möchten und dann diese Kurve im Laufe der Zeit entwickeln sehen.

UPDATE: Least square fitting ist was es ist - ich habe keine bestimmte Routine zu empfehlen. SciPy has one, bin ich mir sicher. Es tut mir leid, dass ich keine bessere Empfehlung habe. Ich lerne jetzt nur Python.

Ich weiß nicht, was Sie mit "Fitness Indikator" meinen. Die Anpassung der kleinsten Quadrate berechnet Koeffizienten, die das mittlere Quadrat des Fehlers zwischen der Anpassung und den Daten an jedem Punkt minimieren.

Nur eine Möglichkeit, mehrere Datensätze zu einer einzigen Anpassung zu kombinieren: Zusammenführen und die Berechnung erneut ausführen.

Answer 3

ich eine umfassende Matlab Toolbox entwickelt passen Parameter und Anfangswerte von Oden an mehreren experimentellen Datensätzen. Je nach Experiment kann es unterschiedliche Anfangswerte verarbeiten und ist unter www.potterswheel.de verfügbar.