In [1]:
import scipy.integrate
import scipy.sparse
import numpy as np
import collections
import matplotlib.pyplot as plt
%matplotlib inline
from copy import deepcopy
from pandas import DataFrame
class IVPResult:
pass
def solve_ivp(f, ts, x0, p=None, integrator='dopri5', store_trajectory=False,
sens=False, dfx=None, dfp=None):
"""
Solve initial value problem
d/dt x = f(t, x, p); x(t0) = x0.
Evaluate Solution at time points specified in ts, with ts[0] = t0.
Compute sensitivity matrices evaluated at time points in ts if sens=True.
"""
if sens:
def f_vode(t, xGxGp, p):
x = xGxGp[:x0.shape[0]]
Gx = xGxGp[x0.shape[0]:x0.shape[0]+x0.shape[0]*x0.shape[0]]
Gx = Gx.reshape([x0.shape[0], x0.shape[0]])
Gp = xGxGp[x0.shape[0]+x0.shape[0]*x0.shape[0]:]
Gp = Gp.reshape([x0.shape[0], p.shape[0]])
dx = f(t, x, p)
dfxev = dfx(t, x, p)
dfpev = dfp(t, x, p)
dGx = dfxev.dot(Gx)
dGp = dfxev.dot(Gp) + dfpev
return np.concatenate(
[dx,
dGx.reshape(-1),
dGp.reshape(-1)])
ivp = scipy.integrate.ode(f_vode)
else:
ivp = scipy.integrate.ode(f)
ivp.set_integrator(integrator)
if store_trajectory:
times = []
points = []
def solout(t, x):
if len(times) == 0 or t != times[-1]:
times.append(t)
points.append(np.copy(x[:x0.shape[0]]))
ivp.set_solout(solout)
if sens:
ivp.set_initial_value(np.concatenate(
[x0,
np.eye(x0.shape[0]).reshape(-1),
np.zeros([x0.shape[0], p.shape[0]]).reshape(-1)
]), ts[0])
else:
ivp.set_initial_value(x0, ts[0])
ivp.set_f_params(p)
result = IVPResult()
result.ts = ts
result.xs = np.zeros([ts.shape[0], x0.shape[0]])
if sens:
result.Gxs = np.zeros([ts.shape[0], x0.shape[0], x0.shape[0]])
result.Gps = np.zeros([ts.shape[0], x0.shape[0], p.shape[0]])
result.success = True
result.xs[0,:] = x0
#meine Vermutung:
if sens:
result.Gxs[0,:,:] = np.eye(x0.shape[0])
for ii in range(1,ts.shape[0]):
ivp.integrate(ts[ii])
result.xs[ii,:] = ivp.y[:x0.shape[0]]
if sens:
result.Gxs[ii,:,:] = ivp.y[x0.shape[0]:x0.shape[0]+x0.shape[0]*x0.shape[0]].reshape([x0.shape[0], x0.shape[0]])
result.Gps[ii,:,:] = ivp.y[x0.shape[0]+x0.shape[0]*x0.shape[0]:].reshape([x0.shape[0], p.shape[0]])
if not ivp.successful():
result.success = False
break
if store_trajectory:
result.trajectory_t = np.array(times)
result.trajectory_x = np.array(points)
return result
In [2]:
def dX_dt(t, X, p): #predator-prey DEs
R = X[0]
F = X[1]
alpha, beta, gamma, delta = p
dRdt = alpha * R - beta * F * R
dFdt = gamma * R * F - delta * F
return np.array([dRdt,dFdt])
def dfx(t, X, p):
R = X[0]
F = X[1]
alpha, beta, gamma, delta = p
return np.array([alpha - beta * F, - beta * R, gamma * F,\
gamma * R - delta]).reshape((2, 2))
def dfp(t, X, p):
R = X[0]
F = X[1]
alpha, beta, gamma, delta = p
return np.array([[R, - F * R, 0., 0.],[0., 0., R * F, -F]])
In [3]:
# delta x should be a vector of dimension equal to the number
# of variables, which is true only if we do not change the parameters.
def evaluate_multiple_shooting(f, dfx, dfp, meas_times, meas_values, s_guesses, p):
F2 = deepcopy(s_guesses[1:])
m, nx = s_guesses.shape
#print((m-1)*nx, m*nx+len(p))
J2 = -1 * np.eye((m-1)*nx, m*nx+len(p), k=nx)
for i in range(m-1):
meas_times_interval = meas_times[i:i+2]
s0 = deepcopy(s_guesses[i])
temp_result = solve_ivp(f, meas_times_interval, s0, p=p, sens=True, dfx=dfx, dfp=dfp)
F2[i] -= temp_result.xs[1]
J2[i*nx:(i+1)*nx,i*nx:(i+1)*nx] = temp_result.Gxs[1]
J2[i*nx:(i+1)*nx,-len(p):] = temp_result.Gps[1]
F1 = meas_values - s_guesses
F1[-1] = meas_values[-1] - temp_result.xs[1]
J1 = np.eye(m * nx, nx*m + len(p))
J1[-nx:,-len(p):] = temp_result.Gps[1]
J1[-nx:,-len(p) - nx: -len(p)] = temp_result.Gxs[1]
return F1.flatten(), -F2.flatten(), -J1, J2
In [4]:
def gauss_newton(F_J, s_guesses, param, meas_times, meas_values, itmax=100, tol=1e-7, verbose=1):
m, nx = s_guesses.shape
tau = 1.
if verbose == 1:
print("Start gauss_newton with itmax=",itmax,"and tol=",tol,"and tau=",tau)
Fs = np.zeros(itmax+1)
i = 0
norm = tol+1.0
dx_prev = 1
while i<itmax and norm > tol:
F1, F2, J1, J2 = F_J(dX_dt, dfx, dfp, meas_times, meas_values, s_guesses, param)
print(DataFrame(J1))
sdf
A = np.vstack((np.hstack((J1.transpose().dot(J1),J2.transpose())),np.hstack((J2,np.zeros((J2.shape[0],J2.shape[0]))))))
b = (-1)*np.hstack((J1.transpose().dot(F1),F2))
Dx = np.linalg.lstsq(A, b)[0][:J1.shape[1]]
param += tau * Dx[-len(param):]
s_increase = Dx[:-len(param)].reshape((m, nx))
norm = np.linalg.norm(Dx)
s_guesses += tau * s_increase
i += 1
print('Iteration:', i)
print('|F1|', np.linalg.norm(F1))
print('|F2|', np.linalg.norm(F2))
print('|Dx|', np.linalg.norm(Dx))
print('|k|', np.linalg.norm(Dx) / np.linalg.norm(dx_prev))
dx_prev = deepcopy(Dx)
'''
print('Dx')
print(Dx)
print('s_incr')
print(s_increase)
print('s_guess')
print(s_guesses)
print('dp')
print(Dx[-len(param):])
print(param)
'''
return s_guesses, param
In [5]:
np.random.seed(31)
meas_times = np.arange(21) * 5
noise = 5. * np.random.rand(21, 2)
s0_init = np.array([20., 10.])
param_init = np.array((.2, .01, .001, .1))
result = solve_ivp(dX_dt, meas_times, s0_init, p=param_init) #just to create the measurements
meas_values = result.xs + noise
print("Initial values close to the correct parameters:")
# We start with an alpha variation ftom the true values
alpha = .01
s_guesses = deepcopy(meas_values)
#as given in the exercise we use the measurements at the grid points as guesses for the initial values
param = param_init * (1 + alpha * np.random.rand(len(param_init.shape)))
s0, param = gauss_newton(evaluate_multiple_shooting, s_guesses,np.zeros(4),meas_times,meas_values)
'''
print('relative error')
print((param - param_init) / param_init)
print((s0 - s0_init)/ s0_init, '\n')
print('estimated parameter vs true ones')
print(param)
print(param_init)
'''
In [6]:
for l in range(20):
x0 = s0[l]
ts = np.array([meas_times[l], meas_times[l+1]])
result = solve_ivp(dX_dt, ts, x0, p=param, integrator='dopri5', store_trajectory=True,
sens=True, dfx=dfx, dfp=dfp)
plt.plot(result.trajectory_t, result.trajectory_x)
In [ ]: