【问题标题】:Python: Unable to solve a differential equation using odeint with signum functionPython:无法使用带有符号函数的 odeint 求解微分方程
【发布时间】:2017-10-06 02:45:55
【问题描述】:

我正在尝试解决这个问题:

U 在哪里

这里:

s=c*e(t)+e_dot(t)

e(t)=theta(t)-thetad(t) 

e_dot(t)=theta_dot(t)-thetad_dot(t)

其中 thetad(theta desired)=sin(t)-- 即要跟踪的信号! 和 thetad_dot=cos(t),J=10,c=0.5,eta=0.5

我首先尝试使用 odeint - 它在 t=0.4 后给出错误,即 theta(上述微分方程的解)下降到 0 并此后保持不变。但是,当我尝试将 mxstep 增加到 5000000 时,我可以在 t=4.3s 之前得到一些正确的图表。

我想得到一个纯正弦图。那就是我希望theta(上述微分方程的解)跟踪thetad,即sin(t)。但在 t=4.3sec 后无法准确跟踪 thetad。在这段时间之后,theta(上述微分方程的解)简单地下降到 0 并停留在那里。

这是我的代码-

from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt

c=0.5
eta=0.5
def f(Y,t):
    x=Y[0]
    y=Y[1]
    e=x-np.sin(t)
    de=y-np.cos(t)
    s=c*e+de
    return [y,-c*(de)-np.sin(t)-eta*np.sign(s)]

t=np.linspace(0,10,1000)
y0=[0.5,1.0]
sol=odeint(f,y0,t,mxstep=5000000)
#sol=odeint(f,y0,t)
print(sol)
angle=sol[:,0]
omega=sol[:,1]
error=angle-np.sin(t)
derror=omega-np.cos(t)
plt.subplot(4,2,1)
plt.plot(t,angle)
plt.plot(t,error,'r--')
plt.grid()
plt.subplot(4,2,2)
plt.plot(t,omega)
plt.plot(t,derror,'g--')
plt.grid()
plt.subplot(4,2,3)
plt.plot(angle,omega)
plt.grid()
plt.subplot(4,2,4)
plt.plot(error,derror,'b--')
plt.savefig('signum_graph.eps')
plt.show()

【问题讨论】:

    标签: python controller ode mode sliding


    【解决方案1】:

    odeint 使用的积分器(可能还有您发现已实现的所有积分器)基于您的导数 (f) 在每个步骤中是连续的(并且具有连续导数)的假设。 signum 函数显然打破了这个要求。这会导致错误估计非常高,无论步长多小,进而导致步长过小以及您遇到的问题。

    对此有两种通用解决方案:

    1. 选择您的步长,使符号翻转正好落在一个台阶上。 (这可能非常乏味。)

    2. 用一个非常陡峭的sigmoid 替换signum 函数,这对于大多数应用程序来说应该没问题,甚至更现实。在您的代码中,您可以使用

      sign = lambda x: np.tanh(100*x)
      

      而不是np.sign。这里的因子 100 控制 sigmoid 的陡度。选择它足够高以反映您实际需要的行为。如果您选择过高,这可能会对您的运行时间产生负面影响,因为积分器会将步长调整为不必要的小。

    在你的情况下,设置一个小的绝对容差似乎也有效,但我不会依赖这个:

    sol = odeint(f,y0,t,mxstep=50000,atol=1e-5)
    

    【讨论】:

    • ODE 函数应该具有中等大小的连续导数,直到方法的阶。什么“中等大小”取决于步长算法的内部结构。甚至阶跃函数的一阶导数也是狄拉克δ峰,越高越差...
    • 非常感谢! Wrzlprmft 你为我省去了很多麻烦!我在Integrate.ode 中查看Vode 选项。不过这里的 tanh 确实帮了我很多忙!
    • 另一个变种是def sign(s): s *= 1e3; return s/(1+np.abs(s));。还要尝试移动 sigmoid 函数以使其不对称,sign = lambda x: np.tanh(100*x-5) 或其他小偏移量。只有当这些变化下的解决方案相似时,您才能希望在给定的区间内存在某种广义的解决方案。
    猜你喜欢
    • 2015-03-05
    • 1970-01-01
    • 2018-02-08
    • 2012-02-28
    • 2015-12-13
    • 2017-11-18
    • 2018-11-01
    • 1970-01-01
    • 1970-01-01
    相关资源
    最近更新 更多