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

Question

我正在尝试制作一个简单的清洁计划工具，用于在热交换器网络中何时进行化学清洁。但是，当我正确地找到最佳清洗时间（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()

Answer 1

如果模拟不需要继续超过该条件，则有一种可变时间方法，将导数除以最终时间变量（tf）。这将调整最终时间，类似于为Jennings benchmark problem显示的方法。下面是一个简化的问题：

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处达到终端约束。换热器问题是否可以采用类似的策略？有模拟换热器清洗的方法，但改变可变时间可能是最简单的解决方案