如何解决为什么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)
解决方法
我无法重现您遇到的问题,因此,我将提供一些具体建议,以对应用程序进行故障排除并提高求解器查找解决方案的能力。
帮助求解器找到解决方案
- 避免使用
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
驱动到目标设定点。
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。希望您会考虑将模型或案例研究贡献给开源存储库。
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。