时间戳

首页 > 解决方案 > 为什么 12 秒后 MHE 对我的模型持怀疑态度？

gekko - 为什么 12 秒后 MHE 对我的模型持怀疑态度？

问题描述

我是控制系统的新手，并尝试使用 GEKKO。对于我的模型模拟，geos 正确并在 MHE 和 MPC 上运行它在 t = 15 秒后正确运行，如果 t = 16 秒在 12 秒，它会如附图所示持怀疑态度。这有什么原因吗？

m = GEKKO(remote = False)

g  = m.Const(value = 9.81)
Cc = m.Const(value = 2*10**-5)
D1 = m.Const(value = 0.1016)
D2 = m.Const(value = 0.1016)
h1 = m.Const(value = 100) 
hv = m.Const(value = 1000)
L1 = m.Const(value = 500)
L2 = m.Const(value = 1100)
V1 = m.Const(value = 4.054)
V2 = m.Const(value = 9.729)
A1 = m.Const(0.008107)
A2 = m.Const(0.008107)

beta1 = m.Const(value = 1.5 * 10**9)
beta2 = m.Const(value = 1.5* 10**9)
Pres = m.Const(value = 1.261 * 10**7)
M = m.Const(value = 1.992 * 10**8) 
rho = m.Const(value = 932)
PI = m.Const(value = 2.32 * 10**(-9))
visc = m.Const(value = 0.024)

fo = m.Const(value = 60)
Inp = m.Const(value = 65)
Pnp = m.Const(value = 1.625 * 10**5)

z = m.MV(1, lb = 0.1, ub = 1)

f = m.MV(35, lb = 35, ub = 65)
f.STATUS = 1
f.DMAX = 30   
f.DCOST = 0.0002 # slow down change of frequency
f.UPPER = 65

Ppin = m.CV()
Ppin.STATUS = 1 # add the SP to the objective
Ppin.SP = 6e6 # set point
Ppin.TR_INIT = 1 # set point trajectory
Ppin.TAU = 5 # time constant of trajectory

Pwf = m.Var(ub = Pres)
Ppout = m.Var()
Pwh = m.Var()
Pm = m.Var()

P_fric_drop = m.Var()
F1 = m.Var()
F2 = m.Var()

q = m.Var(lb = 0)
qc = m.Var(lb = 0)
qr = m.Var(lb = 0)

I = m.Var()
H = m.Var()
BHP = m.Var(lb = 0)
BHPo = m.Var()
qo = m.Var()
Ho = m.Var()

m.Equation([Pwh.dt() == beta2*(q-qc)/V2,
            Pwf.dt() == beta1*(qr-q)/V1,
            q.dt() == (1/M)*(Pwf - Pwh - rho*g*hv - P_fric_drop + (Ppout-Ppin))])

m.Equation([qr == PI*(Pres - Pwf),
            qc == Cc*z*(m.sqrt((Pwh - Pm))),
            Pm == Pwh/2,
            P_fric_drop == F1 + F2,
            F1 == 0.158 * rho * L1* q**2 * visc**0.5 /(D1 * A1**2 * (rho*D1*q)**0.5),
            F2 == 0.158 * rho * L2* q**2 * visc**0.5 /(D2 * A2**2 * (rho*D2*q)**0.5)])

m.Equation([
              qo == q * fo/f,
              H == Ho * (f/fo)**2,
              BHPo == -2.3599e9*qo**3 - 1.8082e7*qo**2 + 4.3346e6*qo + 9.4355e4,
              BHP == BHPo * (f/fo)**3,
              Ppin == Pwf - rho*g*h1 - F1,
              Ho == 9.5970e2+7.4959e3*qo+-1.2454e6*qo**2,
              I == Inp * BHP / Pnp,
              Ppout == H*rho*g + Ppin
             ])

m.options.solver = 1
m.options.MAX_ITER = 250   

m.options.IMODE = 1
m.solve(debug=True)

tf = 30  # final time (sec)
st = 2.0   # simulation time intervals
nt = int(tf/st)+1 # simulation time points
m.time = np.linspace(0,tf,nt)

m.options.CV_TYPE = 2 # squared error
m.options.solver = 3
m.options.NODES=2
m.options.IMODE = 5
m.solve(debug = False)

仿真结果

标签： gekko

解决方案

我无法重现您遇到的问题，因此我将针对您的应用程序故障排除和提高求解器找到解决方案的能力提供一些具体建议。

帮助求解器找到解决方案

  • 避免sqrt(-negative)
#qc == Cc*z*(m.sqrt((Pwh - Pm)))
qc**2 == Cc**2 * z**2 * (Pwh - Pm) # avoid negative in sqrt
  • 避免被零除的可能性。的值q被限制为远离零，但寻找可以通过重新排列来改进的其他方程。
#qo == q * fo/f
qo * f == q * fo

#H == Ho * (f/fo)**2
H * fo**2 == Ho * f**2

#I == Inp * BHP / Pnp
I * Pnp == Inp * BHP
  • 约束 FV 和 MV

它通常有助于缩小DMAX或限制 和 的上限和FVs下限MVs。对Var类型的约束可能导致不可行的解决方案。

  • 降低自由度

对于 MHE 应用程序，我建议您使用FVs而不是MVs用于未知参数。FVs计算所有测量值的单个值。为每个MVs时间点计算一个新的参数值。使用 anFV而不是 anMV可以帮助求解器，因为优化器的决策更少，并且不会调整参数值以跟踪信号噪声。

检查解决方案

  • 配置目标函数。

对于 MHE 应用程序，您需要有一些CV试图匹配的值具有FSTATUS=1. 在您发布的示例中，没有任何插入CVs的FSTATUS=1测量值或任何插入的测量值。CVsMHE 是通过将未知参数调整为FVs或来使模型与测量值对齐MVs。

  • 检查求解器状态以确保它找到了成功的解决方案APPSTAUTS==1
if m.options.APPSTATUS==1:
    print('Successful Solution)
else:
    print('Solver Failed to Find a Solution')
  • 创建变量图，尤其是在有问题的时期。qr您经常可以看到哪里有一个积分变量，例如qc接近一个变量约束。这是我将您的 MHE 应用程序转换为 MPC 应用程序的示例。我正在调整摩擦系数f（不现实）和阀门开z度，以显示 MPC 如何驱动Ppin到目标设定点。

MPC 结果

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt

m = GEKKO(remote = False)

g  = m.Const(value = 9.81)
Cc = m.Const(value = 2*10**-5)
D1 = m.Const(value = 0.1016)
D2 = m.Const(value = 0.1016)
h1 = m.Const(value = 100) 
hv = m.Const(value = 1000)
L1 = m.Const(value = 500)
L2 = m.Const(value = 1100)
V1 = m.Const(value = 4.054)
V2 = m.Const(value = 9.729)
A1 = m.Const(0.008107)
A2 = m.Const(0.008107)

beta1 = m.Const(value = 1.5 * 10**9)
beta2 = m.Const(value = 1.5* 10**9)
Pres = m.Const(value = 1.261 * 10**7)
M = m.Const(value = 1.992 * 10**8) 
rho = m.Const(value = 932)
PI = m.Const(value = 2.32 * 10**(-9))
visc = m.Const(value = 0.024)

fo = m.Const(value = 60)
Inp = m.Const(value = 65)
Pnp = m.Const(value = 1.625 * 10**5)

z = m.MV(1, lb = 0.1, ub = 1)
z.STATUS = 1
z.DMAX = 0.01

f = m.MV(35, lb = 35, ub = 65)
f.STATUS = 1
f.DMAX = 30   
f.DCOST = 0.0002 # slow down change of frequency
f.UPPER = 65

Ppin = m.CV()
Ppin.STATUS = 1 # add the SP to the objective
Ppin.SP = 6e6 # set point
Ppin.TR_INIT = 1 # set point trajectory
Ppin.TAU = 5 # time constant of trajectory

Pwf = m.Var(ub = Pres)
Ppout = m.Var()
Pwh = m.Var()
Pm = m.Var()

P_fric_drop = m.Var()
F1 = m.Var()
F2 = m.Var()

q = m.Var(lb = 0)
qc = m.Var(lb = 0)
qr = m.Var(lb = 0)

I = m.Var()
H = m.Var()
BHP = m.Var(lb = 0)
BHPo = m.Var()
qo = m.Var()
Ho = m.Var()

m.Equation([Pwh.dt() == beta2*(q-qc)/V2,
            Pwf.dt() == beta1*(qr-q)/V1,
            M * q.dt() == Pwf - Pwh - rho*g*hv - P_fric_drop + (Ppout-Ppin)])

m.Equation([qr == PI*(Pres - Pwf),
            qc == Cc*z*(m.sqrt((Pwh - Pm))),
            Pm == Pwh/2,
            P_fric_drop == F1 + F2,
            F1 == 0.158 * rho * L1* q**2 * visc**0.5 /(D1 * A1**2 * (rho*D1*q)**0.5),
            F2 == 0.158 * rho * L2* q**2 * visc**0.5 /(D2 * A2**2 * (rho*D2*q)**0.5)])

m.Equation([
              qo * f == q * fo,
              H * fo**2 == Ho * f**2,
              BHPo == -2.3599e9*qo**3 - 1.8082e7*qo**2 + 4.3346e6*qo + 9.4355e4,
              BHP == BHPo * (f/fo)**3,
              Ppin == Pwf - rho*g*h1 - F1,
              Ho == 9.5970e2+7.4959e3*qo+-1.2454e6*qo**2,
              I * Pnp == Inp * BHP,
              Ppout == H*rho*g + Ppin
             ])

m.options.solver = 1
m.options.MAX_ITER = 250   

m.options.IMODE = 1
m.solve()

tf = 30  # final time (sec)
st = 2.0   # simulation time intervals
nt = int(tf/st)+1 # simulation time points
m.time = np.linspace(0,tf,nt)

m.options.CV_TYPE = 2 # squared error
m.options.solver = 3
m.options.NODES=2
m.options.IMODE = 6

plt.figure(figsize=(10,5))
for i in range(10):
    m.solve(disp=False)
    print(i,m.options.APPSTATUS)
    plt.subplot(10,1,i+1)
    plt.plot(m.time+i*st,Ppin.value)
    plt.xlim([0,50])
    plt.ylim([6e6,7e6])
plt.show()

您的应用程序看起来像是用于钻井中的液压控制。GitHub 上还有其他可用模型。我希望您考虑将您的模型或案例研究贡献给开源存储库。

推荐阅读