如何在GEKKO中实现状态空间模型并完成系泊参数优化

原代码核心问题 #

• GEKKO不支持对数组直接调用.dt()方法，需为每个变量单独定义微分关系
• 错误使用.VALUE提取优化变量初始值修改刚度/阻尼矩阵，应直接将变量嵌入矩阵元素以保留符号关系
• 混用numpy数值函数（如np.cos）与GEKKO符号变量，需替换为GEKKO内置符号函数
• 外力函数需用GEKKO时间参数化方式定义，而非numpy数组
• 约束中使用round()会破坏问题连续性，应直接使用等式约束或设置合理容差

修正后的完整代码 #

import numpy as np
from gekko import GEKKO

# 初始化GEKKO模型
gek = GEKKO(remote=False)

# 仿真参数
max_time_sim = 1600
delta_t = 0.01
t_n = int(max_time_sim / delta_t)
t = gek.Param(value=np.linspace(0, max_time_sim, t_n))

# 物理常数与给定参数
g = 9.81
t_steady = 800

# 假设Tension、stiff、mass_def、damp_def、force_def已提前定义
H = Tension[i-1,0]
T = Tension[i-1,1]
T_max_1 = Tension[i-1,2]  # kN
T_max_23 = Tension[i-1,3] # kN
T_max = T_max_1 + 2 * T_max_23

t_mass = mass_def(T)
c_damp = damp_def(T)
ampphase = force_def(T)
r = 1.0  # 替换为实际给定的r值
r23 = 1.0 # 替换为实际给定的r23值

# 定义优化变量
kx1  = gek.Var(value=1.0, lb=0.0)
cx1  = gek.Var(value=1.0, lb=0.0)
kz1  = gek.Var(value=1.0, lb=0.0)
cz1  = gek.Var(value=1.0, lb=0.0)
kx23 = gek.Var(value=1.0, lb=0.0)
cx23 = gek.Var(value=1.0, lb=0.0)
kz23 = gek.Var(value=1.0, lb=0.0)
cz23 = gek.Var(value=1.0, lb=0.0)

# 构建刚度矩阵（嵌入优化变量）
k_stiff = np.asarray(stiff).copy()
k_stiff[0,0] = gek.Const(k_stiff[0,0]) + kx1 + 2 * kx23
k_stiff[2,2] = gek.Const(k_stiff[2,2]) + kz1 + 2 * kz23
# 将矩阵元素转换为GEKKO常数/变量
for i in range(6):
    for j in range(6):
        if not isinstance(k_stiff[i,j], GEKKO.Var) and not isinstance(k_stiff[i,j], GEKKO.Const):
            k_stiff[i,j] = gek.Const(k_stiff[i,j])

# 构建阻尼矩阵（嵌入优化变量）
damp = np.asarray(c_damp).copy()
damp[0,0] = gek.Const(damp[0,0]) + cx1 + 2 * cx23
damp[2,2] = gek.Const(damp[2,2]) + cz1 + 2 * cz23
# 将矩阵元素转换为GEKKO常数/变量
for i in range(6):
    for j in range(6):
        if not isinstance(damp[i,j], GEKKO.Var) and not isinstance(damp[i,j], GEKKO.Const):
            damp[i,j] = gek.Const(damp[i,j])

# 定义时间相关外力
F = []
for i in range(6):
    amp = gek.Const(H/2 * ampphase[i,0])
    phase = gek.Const(ampphase[i,1])
    F.append(amp * gek.cos(2 * np.pi / T * t + phase))

# 定义运动变量：mot=位移，mot2=速度
mot = gek.Array(gek.Var, (6,1), value=0.0)
mot2 = gek.Array(gek.Var, (6,1), value=0.0)

x    = mot[0,0]
z    = mot[2,0]
rty  = mot[4,0]
dx   = mot2[0,0]
dz   = mot2[2,0]
drty = mot2[4,0]

# 优化约束
gek.Equations([kx1 >= kx23, cx1 >= cx23, kz1 >= kz23, cz1 >= cz23])
gek.Equations([x >= -6.0, x <= 6.0])

# 状态空间方程：速度=位移的微分
gek.Equations([mot2[i,0] == mot[i,0].dt() for i in range(6)])

# 状态空间方程：质量*加速度 + 阻尼*速度 + 刚度*位移 = 外力
for i in range(6):
    acceleration_term = sum(t_mass[i,j] * mot2[j,0].dt() for j in range(6))
    damping_term = sum(damp[i,j] * mot2[j,0] for j in range(6))
    stiffness_term = sum(k_stiff[i,j] * mot[j,0] for j in range(6))
    gek.Equation(acceleration_term + damping_term + stiffness_term == F[i])

# 系缆张力方程
## 系缆1
x2_1 = r * (1 - gek.cos(gek.radians(rty)))
z2_1 = r * gek.sin(gek.radians(rty))
dx2_1 = r * gek.sin(gek.radians(rty)) * gek.radians(drty)
dz2_1 = -r * gek.cos(gek.radians(rty)) * gek.radians(drty)

T1x = kx1 * (x + x2_1) + cx1 * (dx + dx2_1)
T1z = kz1 * (z - z2_1) + cz1 * (dz + dz2_1)
gek.Equation(gek.sqrt(T1x**2 + T1z**2) == T_max_1)

## 系缆2&3
x2_23 = r23 * (1 - gek.cos(gek.radians(rty)))
z2_23 = r23 * gek.sin(gek.radians(rty))
dx2_23 = r23 * gek.sin(gek.radians(rty)) * gek.radians(drty)
dz2_23 = -r23 * gek.cos(gek.radians(rty)) * gek.radians(drty)

T23x = kx23 * (x - x2_23) + cx23 * (dx + dx2_23)
T23z = kz23 * (z + z2_23) + cz23 * (dz + dz2_23)
gek.Equation(gek.sqrt(T23x**2 + T23z**2) == T_max_23)

# 求解设置
gek.options.SOLVER = 1  # 使用APOPT求解器
gek.solve(disp=True)

# 输出优化结果
optimized_params = np.array([kx1.value[0], cx1.value[0], kz1.value[0], cz1.value[0],
                             kx23.value[0], cx23.value[0], kz23.value[0], cz23.value[0]])
print("优化参数：", optimized_params)

关键修正说明 #

1. 状态空间方程实现：通过循环逐个定义变量的微分关系，展开矩阵乘法为逐元素求和，替代原数组级的.dt()调用，匹配12x12矩阵形式的状态空间模型
2. 矩阵变量嵌入：将优化变量直接加入刚度/阻尼矩阵元素，保留符号关系，避免使用.VALUE破坏优化链路
3. 符号函数替换：用gek.cos、gek.radians替代numpy函数，确保GEKKO能进行符号化求解
4. 外力参数化：用GEKKO的Param和cos函数定义时间相关外力，符合GEKKO的符号计算要求
5. 约束修正：移除round()，直接使用等式约束保证张力目标，同时保留位移范围和参数大小约束

内容的提问来源于stack exchange，提问作者Achmad A