我只是试图解决 Schittkowski DAE 测试套件中的 DAE 问题之一(http://klaus-schittkowski.de/mc_dae.htm http://klaus-schittkowski.de/mc_dae.htm)但没有成功(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,100.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)
The resulting plot is
振荡行为将随着计算时间行为的时间步长或搭配节点的增加而减少。另一方面,在 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,0.0,0.0,0.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,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 中似乎没有关于如何解决此类问题的一般建议。我希望有人对此发表评论。
此致,
拉多万