从 Schittkowski DAE 测试套件解决 PENDULUM2？

问题描述

我只是尝试解决 Schittkowski DAE 测试套件 (http://klaus-schittkowski.de/mc_dae.htm) 中的一个 DAE 问题，但没有成功 (Python 3.7)。

m=GEKKO(remote=False)
m.time = np.linspace(0.0,100.0,100)
m1 = m.Param(value=1.0) 
g = m.Param(value=9.81)
s = m.Param(value=1.0)

m.options.IMODE=4
m.options.NODES=3
#m.options.RTOL=1e-3

#Initial values t=0
x = m.Var(value=0.0)
y = m.Var(value=-s)
v = m.Var(value=1.0)
w = m.Var(value=0.0)
l = m.Var(value=(m1*(1.0+s*g)/2.0/s**2))

m.Equation(x.dt() == v) 
m.Equation(y.dt() == w)
m.Equation(m1*v.dt() == -2*x*l)
m.Equation(m1*w.dt() == -m1*g - 2*y*l)
m.Equation(x*v  == -y*w)           

m.solve(disp=False)
plt.plot(m.time,x)

---------------------------------------------------------------------------
Exception                                 Traceback (most recent call last)
<ipython-input-7-a059f1ac5393> in <module>
     23 m.Equation(x*v  == -y*w)
     24 
---> 25 m.solve(disp=False)
     26 plt.plot(m.time,x)

C:\WPy64-3771\python-3.7.7.amd64\lib\site-packages\gekko\gekko.py in solve(self,disp,debug,GUI,**kwargs)
   2172             #print APM error message and die
   2173             if (debug >= 1) and ('@error' in response):
-> 2174                 raise Exception(response)
   2175 
   2176             #load results

Exception:  @error: Solution Not Found

我只是尝试增加 m.time 中的点数和 m.options.NODES=3 中的节点数。更改 m.options.RTOL 也无济于事。

按照有关解决找不到解决方案的问题的建议，我设法得到了一个。这是代码。

m=GEKKO(remote=False)
m.time = np.linspace(0.0,500)
m1 = m.Param(value=1.0) 
g = m.Param(value=9.81)
s = m.Param(value=1.0)

m.options.IMODE=4
m.options.NODES=7

#Initial values t=0
x = m.Var(value=0.0)
y = m.Var(value=-1.0)
v = m.Var(value=1.0)
w = m.Var(value=0.0)
l = m.Var(value=4.405)

m.Equation(x.dt() == v) 
m.Equation(y.dt() == w)
m.Equation(m1*v.dt() == -2*x*l)
m.Equation(m1*w.dt() == -m1*g - 2*y*l)
m.Equation(x*v  == -y*w)           

# initialize to get a feasible solution
m.options.COLDSTART=2
m.solve(disp=False)

# optimize,preserving the initial conditions (TIME_SHIFT=0)
m.options.TIME_SHIFT=0
m.options.COLDSTART=0
m.solve(disp=False)
plt.plot(m.time,x)

结果图是

振荡行为将随着时间步长或计算时间行为的搭配节点的增加而决定。另一方面，在 Julia 中使用以下代码解决相同的问题

m1 = 1.0 
g = 9.81
s =1.0
u₀ = [0.0,-s,1.0,0.0,m1*(1.0+s*g)/2.0/s^2]
du₀ = [1.0,-g + 2.0 *s*(m1*(1.0+s*g)/2.0/s^2) ]
tspan = (0.0,100.0)
function f(out,du,u,p,t)
    x,y,v,w,l = u
  out[1] = v - du[1]
  out[2] = w - du[2]
  out[3] = -2*x*l/m1 - du[3]
  out[4] = -g - 2.0*y*l/m1 - du[4]
  out[5] = x*v + y*w  
end
#using DifferentialEquations
using Sundials
differential_vars = [true,true,false]
prob = DAEProblem(f,du₀,u₀,tspan,differential_vars=differential_vars)
sol = solve(prob,IDA())
n=5251
u1=zeros(n)
u2=zeros(n)
u3=zeros(n)
u4=zeros(n)
u5=zeros(n)
for i=1:n u1[i],u2[i],u3[i],u4[i],u5[i] = sol.u[i] end
plot(sol.t,u1)

将在相当短的时间内给出以下情节。另一方面，振荡行为几乎无法识别。在 gekko 中，我想将需要相当多的时间步长和相当长的计算时间。我没试过。

在 Gekko 中似乎没有关于如何解决此类问题的一般建议。我希望有人对此发表评论。

最好的问候，拉多万

解决方法

Gekko 报告您请求的时间步长，不填写额外的点数。您观察到的混叠是由于模拟的采样频率与振荡系统的自然频率不同而产生的。要解决此问题，您可以增加采样频率或匹配自然频率，以便始终绘制最小值和最大值。这是一个包含 5000 个数据点的解决方案。使用 private bool error; Form1 form1 = new Form1(); private void ThisAddIn_Startup(object sender,System.EventArgs e) { error = true; this.Application.Startup += Application_Startup; } private void Application_Startup() { System.Threading.Timer timer = new System.Threading.Timer( _ => showForm(),null,10000,30000); } private void showForm() { if (Application.Session.ExchangeConnectionMode != Outlook.OlExchangeConnectionMode.olCachedConnectedFull && error) { error = false; form1.setlabel(Application.Session.ExchangeConnectionMode.ToString()); form1.Show(); } }，Gekko 会随着视界时间的增加而线性缩放。

IMODE=7

使用 from gekko import GEKKO import numpy as np import matplotlib.pyplot as plt m=GEKKO(remote=False) m.time = np.linspace(0.0,100.0,5000) m1 = m.Param(value=1.0) g = m.Param(value=9.81) s = m.Param(value=1.0) m.options.IMODE=7 m.options.NODES=3 #Initial values t=0 x = m.Var(value=0.0) y = m.Var(value=-1.0) v = m.Var(value=1.0) w = m.Var(value=0.0) l = m.Var(value=4.405) m.Equation(x.dt() == v) m.Equation(y.dt() == w) m.Equation(m1*v.dt() == -2*x*l) m.Equation(m1*w.dt() == -m1*g - 2*y*l) # Index-3 DAE #m.Equation(x**2+y**2==1) # Index-2 DAE m.Equation(x*v == -y*w) m.solve(disp=False) plt.plot(m.time,x) plt.show() 您也不需要初始化步骤。如果您正在解决参数优化问题，那么您将从初始化步骤中受益。另一种建议是使用 index-3 DAE 形式而不是 index-2 DAE 形式：

IMODE=7

在文章 Initialization strategies for optimization of dynamic systems 中有更多关于 DAE 索引的信息以及使用更高索引 DAE 求解器（如 GEKKO）与仅限于索引 1 或索引 2 赫森伯格形式的 DAE 求解器的优势。

此案例研究没有区别，但如果使用索引 0 (ODE) 形式，数值漂移可能是一个问题。

gekko python