通过同时求解所有时间变化值，在没有时间导数分析表达式的情况下演化动态系统

Question

我有一个动态系统，其中我有一个时间导数残差的表达式，但没有时间导数实际是什么的分析表达式：sum(da_i/dt) ** 2 = f(a) .我可以通过在每个时间步执行优化来解决此问题。我想我可以将其重新转换为最小二乘问题，并在单个求解中的每个时间步长求解每个值。我认为这是一个线性铸造，但我不是 100% 确定。如何形成这个方程组，在一个求解中有效地求解所有时间的值？

'''
my situation is that I know that this relationship is true:
(sum_{i=1}^n d a_i / dt) ** 2 = f(a).
In this example, I assume a is 2 dimensional, f(a) = a1 **2 - a1 * a2,
and I want to evolve the equation for 5 time steps.
'''
def residual(dadt, current_a):
    return (np.sum(dadt) ** 2 - (current_a[0] **2 - np.sum(a) * current_a[0] * current_a[1])) ** 2
'''
My real problem is a little more complicated than this, but the 
basic elements are here. Now, my strategy is to attempt to form 
something like an Ax=b formulation.
'''
'''
Here is an approach that works. I find a values that make the residual zero.
This approach is a little time-consuming because I have to solve the least-squares
problem after each time iteration.
This is a little goofy because I use an optimization engine to solve the problem.
I think it would be better to cast this as a least squares regression problem.
'''
from scipy.optimize import minimize as mini
n_time_steps = 300
a0 = np.ones(2) * .5
a = a0
dt = 1e-3
a0_record = []
for _ in range(n_time_steps):
    a0_record.append(a[0])
    dadt = mini(residual, np.ones(a.size), args=(a,)).x
    a += dadt * dt
# Let's plot the time history. Maybe that can help verify the new approach.
import matplotlib.pyplot as plt
plt.plot(a0_record)
plt.savefig('a0baseline.pdf')
plt.clf()

'''
I want to cast this into a single least-squares problem. In this case, I think
I should have 2 * 300 = 600 free variables. Since I only have 300 time steps
to enforce the residual to be zero, I think that my system is underdetermined.
I don't quite know how to do this. I sketched my best shot so far below but I 
am sure that there is much improvement to make.
'''
a = np.ones((2, n_time_steps))
A = np.stack([residual((a[i] - a[i-1]) / dt, a[i-1]) for i in range(1,n_time_steps)])
xx = np.linalg.lstsq(A, np.zeros(n_time_steps))

Answer 1

这不是

一个线性问题。您所要做的就是使用 least_squares 或其他非线性求解器来解决它，或者使用隐式方法。请注意，此显式方法需要很长时间才能运行。

from scipy.optimize import least_squares
def fitness(A):
    mysum = 0.
    A = A.reshape(n_time_steps, 2)
    for ii in range(n_time_steps-1):
         mysum += residual((A[ii+1] - A[ii]) / dt, A[ii]) ** 2
    return np.array(mysum)
least_squares(fitness, np.random.random(600)) # beware of the trivial solution a=0