通过在同一时刻求解所有的时变值，发展动态系统而不需要时间导数的解析表达式

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