Почему MHE через 12 секунд скептически относится к моей модели?

ProgramBox

Почему MHE через 12 секунд скептически относится к моей модели?

Post author:admin
Запись опубликована:15 марта, 2023
Post category:Вопросы по программированию

#gekko

#гекко

Вопрос:

Я новичок в системах управления и пытался использовать 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)

Ответ №1:

Я не смог воспроизвести проблему, с которой вы столкнулись, поэтому я дам несколько конкретных советов по устранению неполадок в вашем приложении и улучшению способности решателя находить решение.

Помогите решателю найти решение

Избегайте 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 Вычисляет новое значение параметра для каждого момента времени. Использование an FV вместо an 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()

Похоже, что ваше приложение предназначено для управления гидравлическим давлением при бурении. На GitHub доступны дополнительные модели. Я надеюсь, что вы рассмотрите возможность добавления своей модели или тематических исследований в репозиторий с открытым исходным кодом.

Вопрос:

Ответ №1:

Вам также может понравиться

Попытка индексировать данные в Amazon elastic search

Laravel voyager с phpunit

Razor / Syncfusion: при асинхронной загрузке страницы древовидные данные не будут заполняться, если я не выполню действие