Gekko中DAE边值问题控制变量平滑化方案咨询

我正在求解一个涉及三个控制变量的DAE边值问题，通过下方最小工作示例（MWE），需利用这三个控制变量对最终时刻的8个未知变量进行优化。但得到的控制变量存在剧烈波动，不适用于实际控制。尝试通过DCOST选项施加惩罚项后，求解无法收敛，请问还有其他可行方法吗？

from gekko import GEKKO
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import math
from scipy import special

tg=10
sample_rate= 100
tlist= np.linspace(0,tg,int(tg*sample_rate))

m = GEKKO()
m.time= tlist

# 给定时间函数
sigma= 1.25

def eps_G(t):
    A= np.pi/(np.sqrt(2*np.pi)*sigma*special.erf(tg/(2*np.sqrt(2)*sigma))- np.exp(-(0-tg/2)**2/(2*sigma**2))*tg)
    B= A*np.exp(-(0-tg/2)**2/(2*sigma**2))
    return  A*np.exp(-((t-(tg/2))/sigma)** 2/2)-B

omg_0= [eps_G(t) for t in tlist]
omg_1= np.array(omg_0)
omg= m.MV(omg_1)
omg.STATUS=0

# 控制变量
omx= m.MV(value=0, lb= -2, ub= 2)
omy= m.MV(value=0, lb= -2, ub= 2)
det= m.MV(value=0, lb= -2, ub= 2)

omx.STATUS= 1
omy.STATUS= 1
det.STATUS= 1

# omx.DCOST=0.00001
# omy.DCOST=0.00001
# det.DCOST=0.00001

# 构造方程的numpy数组
Hz= np.array([[0,0],[0,1]])
Hx= np.array([[0,1],[1,0]])
Hy= np.array([[0,-1],[1,0]])

HRR= np.empty((2,2), dtype= object)
for i in range(2):
    for j in range(2):
        HRR[i,j]= det*Hz[i,j]+ (omx/2)*Hx[i,j]

HRI= np.empty((2,2), dtype= object)
for i in range(2):
    for j in range(2):
        HRI[i,j]=  (omy/2)*Hy[i,j]

# 未知变量
DR = m.Array(m.Var,(2,2), value=0, lb= -10, ub= 10)
for i in range(2):
    for j in range(2):
        if i==j:
            DR[i,j].value = 1
        else:
            DR[i, j].value= 0

DI = m.Array(m.Var,(2,2),value=0, lb= -10, ub= 10)
for i in range(2):
    for j in range(2):
        DI[i, j].value=0

# 方程
m.Equations([np.transpose(DR)[i,j]*DR[i,j].dt()+ np.transpose(DI)[i,j]*DI[i,j].dt()== np.dot(np.transpose(DR),np.dot(HRR, DI))[i,j]+np.dot(np.transpose(DR),np.dot(HRI, DR))[i,j]- np.dot(np.transpose(DI),np.dot(HRR, DR))[i,j]+np.dot(np.transpose(DI),np.dot(HRI, DI))[i,j] for i in range(2) for j in range(2) ])

m.Equations([np.transpose(DR)[i,j]*DI[i,j].dt()- np.transpose(DI)[i,j]*DR[i,j].dt()== (omg/2)*Hx[i,j]- np.dot(np.transpose(DR),np.dot(HRR, DR))[i,j]+np.dot(np.transpose(DR),np.dot(HRI, DI))[i,j]- np.dot(np.transpose(DI),np.dot(HRR, DI))[i,j]-np.dot(np.transpose(DI),np.dot(HRI, DR))[i,j] for i in range(2) for j in range(2) ])

# 最终时刻参数
p = np.zeros(int(tg*sample_rate))
p[-1] = 1.0
final = m.Param(value=p)

# 目标函数
# m.Obj(omx*final)
# m.Obj(omy*final)
# m.Obj(det*final)

for i in range(2):
    for j in range(2):
        m.Obj((DR[i,j]- np.identity(2)[i,j])*final)
        m.Obj(DI[i,j]*final)

# 求解DAE
m.options.SOLVER= 3
m.options.IMODE = 6
m.options.MAX_ITER=500
m.solve(disp= True, debug=0)

可行解决方法 #

• 限制控制变量变化速率：给每个控制MV变量添加DMAX参数，直接约束单步变化幅度，例如：

  omx.DMAX = 0.1  # 控制变量omx每步变化不超过0.1
  omy.DMAX = 0.1
  det.DMAX = 0.1
  

  这种硬约束比DCOST更直接，能有效抑制剧烈波动。
• 自定义平滑惩罚项：在目标函数中添加控制变量导数的平方项，通过调整权重平衡优化目标和平滑性，避免DCOST的收敛问题：

  # 添加控制变量变化量的平方惩罚
  m.Obj(1e-4 * (omx.dt())**2)
  m.Obj(1e-4 * (omy.dt())**2)
  m.Obj(1e-4 * (det.dt())**2)
  

  可根据实际情况调整1e-4权重值，权重越小平滑约束越弱，反之越强。
• 优化时间网格：
  • 降低采样率（如sample_rate=20），减少问题规模，提升求解器收敛性；
  • 采用非均匀时间网格，在函数变化剧烈的关键时段加密采样，其他时段稀疏采样，例如：

    # 前2秒密采样，后8秒疏采样
    t1 = np.linspace(0,2,40)
    t2 = np.linspace(2,10,20)
    tlist = np.concatenate((t1,t2))
    m.time = tlist
    

• 调整求解器与参数：
  • 切换求解器：将m.options.SOLVER=3（IPOPT）改为m.options.SOLVER=1（APOPT），APOPT在处理非线性DAE和优化问题时稳定性更好；
  • 放宽容忍度：增大m.options.OTOL（优化容忍度）或m.options.RTOL（残差容忍度），例如设置为1e-4，降低求解难度；
• 优化初始轨迹：给控制变量设置平缓的初始轨迹（如线性变化曲线），而非全0初始值，帮助求解器找到更优起始点；
• 约束二阶导数：如果需要更平滑的轨迹，可添加控制变量二阶导数的约束，进一步限制剧烈波动：

  # 限制omx的二阶导数绝对值不超过0.05
  m.Equation(omx.dt().dt() <= 0.05)
  m.Equation(omx.dt().dt() >= -0.05)
  # 同理设置omy和det的二阶导数约束
  

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