Gekko求解非线性方程时'results.json'未找到错误求助

错误本质 #

这个报错并非文件缺失问题，而是求解器因模型异常（方程无解、数值发散、变量/参数定义错误）提前终止，导致无法生成结果文件。结合你修改的t、F2、theta1变量，重点排查以下代码问题：

代码修复要点 #

1. 混淆Param与Var的定义 #

你将依赖动态变量的tau1、tau2、s1等声明为m.Param()（参数，固定值），但实际它们是随变量变化的动态值，必须改为m.Var()并补充方程定义：

# 替换原Param声明为Var
tau1 = m.Var()
tau2 = m.Var()
tau3 = m.Var()
tau4 = m.Var()
s1 = m.Var()
s2 = m.Var()
s3 = m.Var()
s4 = m.Var()

# 补充动态计算方程
m.Equation(tau1 == m.abs(m.atan((w/2-y_p)/(x_p+l/2))))
m.Equation(tau2 == m.abs(m.atan((w/2-y_p)/(l/2-x_p))))
m.Equation(tau3 == m.abs(m.atan((y_p+w/2)/(x_p+q/2))))
m.Equation(tau4 == m.abs(m.atan((y_p+w/2)/(q/2-x_p))))
m.Equation(s1 == m.sqrt((l/2+x_p)**2 + (w/2-y_p)**2))
m.Equation(s2 == m.sqrt((l/2-x_p)**2 + (w/2-y_p)**2))
m.Equation(s3 == m.sqrt((q/2+x_p)**2 + (w/2+y_p)**2))
m.Equation(s4 == m.sqrt((q/2-x_p)**2 + (w/2+y_p)**2))

2. 参数赋值方式错误 #

F、F2、F3作为m.Param()，不能用普通赋值语句绑定，需改为Gekko认可的初始化或方程形式：

# 直接初始化固定参数
F = m.Param(mass*g)
F3 = m.Param(mass*g*0.7)
# 条件表达式的F2需用Var+方程定义
F2 = m.Var()
m.Equation(F2 == m.if2(t-0.125, mass*g*0.7, 0))

3. 重复覆盖Gekko变量 #

theta5先声明为m.Param()，后又用普通赋值覆盖，导致模型逻辑断裂，改为动态变量+方程：

theta5 = m.Var()
m.Equation(theta5 == 2*math.pi/3*t)

4. 数值稳定性优化 #

• 5000个时间步会大幅增加计算量，先减少到500步测试收敛性：

  m.time=np.linspace(0.046875,0.328125,500)
  

• 给非线性变量添加初始值和上下限，帮助求解器收敛：

  theta1 = m.Var(value=-1.89/180*math.pi, lb=-math.pi/2, ub=math.pi/2)
  

5. 开启调试输出 #

关闭disp=False，改为disp=True，查看求解器的详细报错信息，直接定位收敛失败的具体方程：

m.solve(disp=True)

核心修正后代码片段 #

from gekko import GEKKO
import numpy as np
import math
import matplotlib.pyplot as plt
import pandas as pd

m = GEKKO(remote=False)
# 先减少时间步长测试
m.time=np.linspace(0.046875,0.328125,500)

# Variables
xc0 = 0
yc0 = 0
k3_ratio= 1.0
strightness = 0e-6
x_p = m.Var(0)
y_p = m.Var(0.5*(1-k3_ratio)/(1+k3_ratio))
theta = m.Var(value=0)
x_c = m.Var(xc0)
y_c = m.Var(yc0)
a = m.Var(0)
b = m.Var(0)
z = m.Var(0)
ac = m.Var(0)
bc = m.Var(0)
t = m.Param(m.time)
R1 = m.Var()
R2 = m.Var()
R3 = m.Var()
R4 = m.Var()

# 修正为Var并后续绑定方程
tau1 = m.Var()
tau2 = m.Var()
tau3 = m.Var()
tau4 = m.Var()
s1 = m.Var()
s2 = m.Var()
s3 = m.Var()
s4 = m.Var()

theta5 = m.Var()
theta1 = m.Var(value=-1.89/180*math.pi)
theta3 = m.Var()
A = m.Var()
B = m.Var()
C = m.Var()
K = m.Var()

theta11 = m.Var()
theta33 = m.Var()
AA = m.Var()
BB = m.Var()
CC = m.Var()
KK = m.Var()

theta333 = m.Var()
theta3333 = m.Var()
d = m.Var()
d2 = -0.060904
d3 = 31.521

# 修正参数定义
F = m.Param(mass*g)
F2 = m.Var()
F3 = m.Param(mass*g*0.7)

# constants
h = 3.78/180*math.pi
beta = 45/360*2*math.pi
r1 = 379
r2 = 188
r3 = 70
r0 = 486.6
r11 = 140
r22 = 519.299
r33 = 58
r00 = 526.647
l = 0.233
q = l
w = 0.119
mu = 0.02
g = 9.8
k1 = 2.5e8
k2 = k1
k3 = k1 *k3_ratio
k4 = k3
mass = l*w*0.012*2790
theta0 = 22.728/360*2*math.pi
theta00 = 80.418/360*2*math.pi
c1 = 10
c2 = c1
c3 = c1 *k3_ratio
c4 = c3

# 补充tau和s的方程
m.Equation(tau1 == m.abs(m.atan((w/2-y_p)/(x_p+l/2))))
m.Equation(tau2 == m.abs(m.atan((w/2-y_p)/(l/2-x_p))))
m.Equation(tau3 == m.abs(m.atan((y_p+w/2)/(x_p+q/2))))
m.Equation(tau4 == m.abs(m.atan((y_p+w/2)/(q/2-x_p))))
m.Equation(s1 == m.sqrt((l/2+x_p)**2 + (w/2-y_p)**2))
m.Equation(s2 == m.sqrt((l/2-x_p)**2 + (w/2-y_p)**2))
m.Equation(s3 == m.sqrt((q/2+x_p)**2 + (w/2+y_p)**2))
m.Equation(s4 == m.sqrt((q/2-x_p)**2 + (w/2+y_p)**2))

# 剩余方程保持不变
m.Equation(theta5 == 2*math.pi/3*t)
m.Equation(theta1 == -1.89/180*math.pi+h*(0.28005+0.43989*theta5/beta-0.315055*m.cos(4*math.pi/3*theta5/beta-math.pi/6)))
m.Equation(A == 2*r0*r3*m.cos(theta0)-2*r1*r3*m.cos(theta1))
m.Equation(B == 2*r0*r3*m.sin(theta0)-2*r1*r3*m.sin(theta1))
m.Equation(C == r0**2+r1**2+r3**2-r2**2-2*r0*r1*(m.cos(theta0)*m.cos(theta1)+m.sin(theta0)*m.sin(theta1)))
m.Equation(K == (-B-(B**2-C**2+A**2)**(1/2))/(C-A))
m.Equation(theta3 == m.atan(K)*2)
m.Equation(theta11 == theta3 + math.pi)
m.Equation(AA == 2*r00*r33*m.cos(theta00)-2*r11*r33*m.cos(theta11))
m.Equation(BB == 2*r00*r33*m.sin(theta00)-2*r11*r33*m.sin(theta11))
m.Equation(CC == r00**2+r11**2+r33**2-r22**2-2*r00*r11*(m.cos(theta00)*m.cos(theta11)+m.sin(theta00)*m.sin(theta11)))
m.Equation(KK == (-BB+(BB**2-CC**2+AA**2)**(1/2))/(CC-AA))
m.Equation(theta33 == m.atan(KK)*2)
m.Equation(theta333 == theta33.dt())
m.Equation(theta3333 == theta333.dt())
m.Equation(d == 0.145209 - m.cos(theta33)*0.116)
m.Equation(F2 == m.if2(t-0.125, mass*g*0.7, 0))

# Reaction force及动力学方程保持不变
m.Equations([z == theta.dt(), a== x_p.dt(), b== y_p.dt()])
m.Equations([ac== x_c.dt(), bc== y_c.dt()])
# 此处省略原R1-R4的长方程，保持原内容不变
m.Equation(x_c == yc0/(w/2)*(-s1*m.cos(tau1+theta)+s2*m.cos(tau2-theta)+s3*m.cos(tau3-theta)-s4*m.cos(tau4+theta))/4 +(4*x_p-s1*m.cos(tau1+theta)+s2*m.cos(tau2-theta)-s3*m.cos(tau3-theta)+s4*m.cos(tau4+theta) )/4)
m.Equation(y_c == yc0/(w/2)*( s1*m.sin(tau1+theta)+s2*m.sin(tau2-theta)+s3*m.sin(tau3-theta)+s4*m.sin(tau4+theta) )/4 +(4*y_p+s1*m.sin(tau1+theta)+s2*m.sin(tau2-theta)-s3*m.sin(tau3-theta)-s4*m.sin(tau4+theta))/4 )
m.Equation(F3-mu*mass*g-mu*(R1+R2+R3+R4+k4*m.abs(strightness))==mass*a.dt())
m.Equation(-R1+R2-R3+R4-mu*mass*g-k4*strightness-F-F2==mass*b.dt())
m.Equation(F3*(d-y_p)-(x_p-mu*y_p+(l+mu*w)/2)*R1 +(x_p+mu*y_p-(l+mu*w)/2)*R2 -(x_p-mu*y_p+(q-mu*w)/2)*R3 + (x_p+mu*y_p-(q-mu*w)/2)*R4+k4*strightness + F*(d2-x_p) + F2*(d3-x_p)== mass*((w**2+l**2)/12+(x_p-x_c)**2+(y_p-y_c)**2)*z.dt())

m.options.IMODE=6
m.solve(disp=True)
df = pd.DataFrame([m.time, theta.value])
df.to_csv('D:/pyyhon/1축 힘 변화X/yaw angle_middle.csv',index=None)
plt.plot(m.time, theta.value)
plt.show()

内容的提问来源于stack exchange，提问作者Ginie Kim