为什么12秒钟后MHE怀疑地运行在我的模型上？

编程问答 2022-06-11

问题描述

我是控制系统的新手，并尝试使用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)

Simulation results

解决方法

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

帮助求解器找到解决方案

  • 避免使用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为每个时间点计算一个新的参数值。使用FV代替MV可以帮助求解器，因为优化器的决策较少，并且不会调整参数值来跟踪信号噪声。

检查解决方案

  • 配置目标函数。

对于MHE应用程序，您需要具有匹配的CV的一些FSTATUS=1值。在您发布的示例中，没有CVs和FSTATUS=1或任何插入的度量。 MHE可以通过将未知参数调整为CVs或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 Results 

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()

您的应用程序看起来像是用于钻井中的液压控制。有additional models available on GitHub。希望您会考虑将模型或案例研究贡献给开源存储库。

gekko