满足语句时将变量设置为零的问题

Posted 2023-02-18

【中文标题】满足语句时将变量设置为零的问题

【英文标题】：Problem setting variable to zero when statement is met

【发布时间】：2021-11-17 16:59:48

【问题描述】：

我正在尝试制作一个简单的清洁计划工具，用于确定何时在热交换器网络中进行化学清洁。但是当我正确地找到清洁的最佳时间（x 变量）时，我无法在时间 t 将缩放厚度设置为零（sigma），我尝试使用 m.if3 但无济于事。我在下面添加了我的问题的简单版本。感谢您提供任何反馈。

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
def LN(x):
            return m.log(x)/np.log(2.718)
m = GEKKO(remote=False)
lambdag=0.1 #[W/mK]
days_to_consider = 1
m.time=np.linspace(0, 24*days_to_consider, 24*days_to_consider+1)
N = 6 #Number of heat exchanger
sigm = m.Array(m.Var,N,value=0.0,lb=0)
Rf = m.Array(m.Var,N,value=0.0,lb=0) #[m2K/W]
U = m.Array(m.Param,N,lb=0)
LMTD = m.Array(m.Param,N,lb=0)
Tco = m.Array(m.Param,N,lb=0)
Tci = m.Array(m.Param,N,lb=0)
Q = m.Array(m.Param,N,value=0.0)
dQ = m.Array(m.Var,N,value=0.0)
x = m.Array(m.MV,N,value=0,lb=0,ub=1,integer=True)
x[0].STATUS=1
x[1].STATUS=1
x[2].STATUS=1
x[3].STATUS=1
x[4].STATUS=1
x[5].STATUS=1
EL = m.Array(m.Param,N,value=0)
ELchc = m.Array(m.Param,N,value=0)

Thilist = [105,116,125,129,136,142] #Hot vapor entering [degC] ->Condensing
mdotlist = [582.5,582.5,582.5,582.5,582.5,582.5] # Solution flow [t/h]
Arealist = [600,400,200,300,200,300] #Heating surface [m2]
kglist = [0.0094,0.0003,0.0007,4.5019e-05,0.0003,4.6977e-05] # Deposit rate
Ucllist = [1700,2040,3300,3300,3200,2300] # Cleaned Heat transfer Coefficient [W/m2K]
Qcllist = [10036.4,9336.6,7185.8,5255.4,5112.5,5678.8]
CE = 0.5 #fuel cost[EUR/kWh]
Cchc = 500 #Cleaning cost [EUR/CIP]
#Temperature into heat exchanger network (HEN)
Tci[0] = 90 # degC
#Loop through HEN
for u in range(0,N):
    Thi = Thilist[u]
    Tci = Thi-8 
    mdot = mdotlist[u]
    Area=Arealist[u]
    # Scaling kinematics
    kg = kglist[u]
    Ucl = Ucllist[u]
    Qcl = Qcllist[u]
    m.Equation(sigm[u].dt()==kg*lambdag)
    #TODO PROBLEM: cannot set sigma to zero at time t when x(t) is 1
    #b = m.if3(x[u]-1,1,0)         # binary switch
    m.Equation(sigm[u]==(1)*Rf[u]*lambdag) 
    U[u] = m.Intermediate(Ucl/(1+Ucl*Rf[u]))
    # Thermodynamics 
    LMTD[u]=m.Intermediate(((Thi-Tci)-(Thi-Tco[u]))/LN((Thi-Tci)/(Thi-Tco[u])))
    Tco[u]=m.Intermediate(LMTD[u]*U[u]*Area/(mdot/3.6*3300*1000)+Tci)
    Q[u]=m.Intermediate(U[u]*Area*LMTD[u]/1000)  
    m.Equation(dQ[u].dt()==1/6*(Qcl - Q[u])) 
    EL[u]=m.Intermediate(CE*dQ[u]) 
    ELchc[u]=m.Intermediate(CE*(Q[u] -1/6*Q[u] )*2.44+Cchc) 
    u +=1 
m.Minimize(m.sum([EL[u]*(1-x[u])+(ELchc[u]*x[u]) for u in range(0,len(x))]))
#Constrains
m.Equation(m.sum(x)<=1.0) # Only one clean at time

m.options.IMODE=6
m.solver_options = ['minlp_maximum_iterations 500', \
                                            'minlp_gap_tol 0.01',\
                                            'nlp_maximum_iterations 500']
m.options.SOLVER = 1
m.solve(debug=True,disp=True)

plt.figure(figsize=(12, 6))
plt.subplot(141)
for i in range(0,5):
    plt.bar(m.time,x[i].value,label='CIP'+str(i), width=1.0)
plt.legend()
plt.subplot(142)
plt.plot(m.time,EL[0].value,label='Energy cost')
plt.plot(m.time,ELchc[0].value,label='CIP cost')
plt.legend()
plt.subplot(143)
for i in range(0,5):
    plt.plot(m.time,U[i].value,label='U'+str(i))
plt.legend()
plt.subplot(144)
for i in range(0,5):
    plt.plot(m.time,sigm[i].value,label='scaling'+str(i))
plt.legend()
plt.show()

【参考方案1】：

如果模拟不需要继续超过该条件，则可以使用可变时间方法将导数除以最终时间变量 (tf)。这将调整最终时间，类似于 Jennings benchmark problem 的方法 show。这是一个简化的问题：

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

m = GEKKO()

m.time = np.linspace(0,1,101)
x = m.Var(0,lb=0,ub=1)
tf = m.FV(1,lb=0.1,ub=10); tf.STATUS=1
m.Equation(x.dt()/tf==0.2)
m.Maximize(tf)

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

t = [ti*tf.value[0] for ti in m.time]

plt.plot(t,x.value,label='x')
plt.legend()
plt.show()

最后时间调整为最大化tf，同时遵守x<1的约束。微分方程为dx/dt=0.2，因此在t=5 处达到了终端约束。可以对热交换器问题使用类似的策略吗？有一些方法可以模拟热交换器的清洗，但更改可变时间可能是最简单的解决方案。

【讨论】：

我试图模拟的是一个超过 100 天的清洁周期，其中决策变量 x 是一个从 0 到 1 的整数，1 是执行清洁过程从而去除所有污垢/结垢 I因此需要将缩放厚度“重置”为 0。问题是 MINLP 问题。

没有办法任意设置sigm[u]的初始条件，因为微分方程m.Equation(sigm[u].dt()==kg*lambdag决定了解。您可以尝试在污染方程上使用“清洁”术语，例如：m.Equation(sigm[u].dt()==kg*lambdag-clean_rate*b)，当满足条件时减少污染。

是的，这也是我得出的结论 :) 约翰，你知道如何确保 xi=1 的时间段不能超过让我们说 2 小时吗？我认为您的第一个答案可能是一个解决方案，但很难制定有用的东西:(

可以使用 m.Equation(m.integral(xi/dt)<2) 之类的东西强制执行基于时间的条件约束。