使用Python模拟批处理操作的推荐方法是什么，特别是在恢复过程中需要将累加器重置为0的情况下。

Question

我们正在探索n列列车的生产调度系统。每列列车在再生几个小时之前，最大吞吐量可达5800次。当在上时，可以将可变的生产速率应用于每一列列车。理想情况下，一次只能有一列火车处于再生状态。目标是按照一个目标函数来协调生产速度和再生时间，通常是地平线上的总生产量。

当我努力为5列列车获得一个可行的解决方案时，我将问题简化为一列列车，其临时目标是测试再生调度。这个问题目前是离散的，因此整个再生周期适合于一个时间步骤。代码附在下面。

# -*- coding: utf-8 -*-
"""
Created on Tue Aug 23 20:37:26 2022

@author: Jacques

Single Train simplification of a more general problem of multi-train batch operation
Train is to be operated in batch mode.  Once total accumulated throughput of
max 5800 is achieved, it enters regeneration that effectively resets accumulated
throughput to 0 after which the cycle starts again. 

The single train case is use to test certain modelling constructs before the problem
is more generalised to more parallel trains and other constraints.

"""

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

m=GEKKO(False)                  #instantiate GEKKO object
m.time=np.linspace(0,23,24)     #define timespace

#Define production flow rates, when on, for each train
Train1_F=m.Var(value=300,lb=100,ub=350,name='Train1F')  #Train 1 Flow

#Helper variable representing regeneration rate
RegenRate1=m.Var(lb=0,ub=5800,name='RegenRate1')  #RegenRate1

#define regen state as binary variable --  1 is online, 0 is regeneration
Train1_On=m.Var(value=1,lb=0,ub=1,integer=True,name='Train1 Online')

#define production variable as product of regen status and production rate - i.e 0 when on regen
Train1_Prod=m.Intermediate(Train1_F*Train1_On,name='Train1 Prod')

#maintain production accumulation per train,  variable per train,random starting points
Train1_Acc=m.Var(value=2000,lb=0,ub=5800,name='Train1_Acc')


#define accumulation as function of production - reset when regen is on
m.Equation(Train1_Acc.dt()==4*Train1_Prod - RegenRate1*(1-Train1_On))

#maximise production
Total_production=m.Intermediate(Train1_Prod)

#Attempts to Maximise production with as little as possible regen-cycles
m.Maximize(Total_production-0.1*(1-Train1_On))

m.options.IMODE=6  #control mode
m.options.SOLVER=1  #select IPOT, suitable for Mixed Integer Programming
m.options.NODES=2
m.solve()

plt.plot(Train1_Acc)
plt.show()
plt.plot(Train1_On)

我想提请注意Train1_Acc的定义，在这个定义中，我试图使积累成为一个连续的函数，即在列车为OFF时，使积累率与生产率成正比，并且具有足够大的任意值，以表示重新设置吞吐量积累的再生步骤。

求解后，该列车的吞吐量图如下：

​

​

值得注意的是，最后两个循环的再生不是0--大概是因为它在目标函数的最终值方面是不重要的。实际上，在我们的背景下，这样的解决方案将是模棱两可的或错误的解释。我们如何用模型的重置部分重新表述累积，在每个再生步骤中强制显式重置为零？

Answer 1

使用MV (操纵变量)作为决策Train1_On而不是Var。这允许变量有一个更改成本DCOST，这将惩罚再生周期不必要的移动。

#define regen state as binary variable --  1 is online, 0 is regeneration
Train1_On=m.MV(value=1,lb=0,ub=1,integer=True,name='Train1 Online')
Train1_On.STATUS=1
Train1_On.DCOST =1.0

另一个建议是使用最大的再生率和最小的再生率时，不再生。这提高了求解速度，因为目标函数驱动该值，而不是让这些变量不确定。

m.Maximize(10*Total_production)
m.Maximize(100*RegenRate1_Off)
m.Minimize(RegenRate1_On)
m.Minimize(0.1*(1-Train1_On))

​

​

解决方案要快得多，并且在大多数循环中都达到零。附加的目标函数权重可能会有一些改进。

# -*- coding: utf-8 -*-
"""
Created on Tue Aug 23 20:37:26 2022

@author: Jacques

Single Train simplification of a more general problem of multi-train batch operation
Train is to be operated in batch mode.  Once total accumulated throughput of
max 5800 is achieved, it enters regeneration that effectively resets accumulated
throughput to 0 after which the cycle starts again. 

The single train case is use to test certain modelling constructs before the problem
is more generalised to more parallel trains and other constraints.

"""

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

m=GEKKO(remote=False)           #instantiate GEKKO object
m.time=np.linspace(0,23,24)     #define timespace

#Define production flow rates, when on, for each train
Train1_F=m.Var(value=300,lb=100,ub=350,name='Train1F')  #Train 1 Flow

#Helper variable representing regeneration rate
RegenRate1=m.Var(lb=0,ub=5800,name='RegenRate1')  #RegenRate1

#define regen state as binary variable --  1 is online, 0 is regeneration
Train1_On=m.MV(value=1,lb=0,ub=1,integer=True,name='Train1 Online')
Train1_On.STATUS=1
Train1_On.DCOST =1.0

#define production variable as product of regen status and production rate - i.e 0 when on regen
Train1_Prod=m.Intermediate(Train1_F*Train1_On,name='Train1 Prod')

#maintain production accumulation per train,  variable per train,random starting points
Train1_Acc=m.Var(value=2000,lb=0,ub=5800,name='Train1_Acc')
RegenRate1_Off = m.Intermediate(RegenRate1*(1-Train1_On))
RegenRate1_On = m.Intermediate(RegenRate1*Train1_On)

#define accumulation as function of production - reset when regen is on
m.Equation(Train1_Acc.dt()==4*Train1_Prod - RegenRate1_Off)

#maximise production
Total_production=m.Intermediate(Train1_Prod)

#Attempts to Maximise production with as little as possible regen-cycles
m.Maximize(10*Total_production)
m.Maximize(100*RegenRate1_Off)
m.Minimize(RegenRate1_On)
m.Minimize(0.1*(1-Train1_On))

m.options.IMODE=6   #control mode
m.options.SOLVER=1  #select APOPT, suitable for Mixed Integer Programming
m.options.NODES=2
m.solve()

plt.figure(figsize=(12,5))
plt.subplot(2,1,1)
plt.plot(m.time,Train1_Acc,'b.-',label='Acc')
plt.plot(m.time,RegenRate1,'g.-',label='Regen')
plt.legend(); plt.grid()
plt.subplot(2,1,2)
plt.step(m.time,Train1_On,'r.-',label='Train On')
plt.legend(); plt.grid(); plt.xlabel('Time')
plt.show()

 Number of state variables:    207
 Number of total equations: -  138
 Number of slack variables: -  0
 ---------------------------------------
 Degrees of freedom       :    69
 
 ----------------------------------------------
 Dynamic Control with APOPT Solver
 ----------------------------------------------
Iter:     1 I:  0 Tm:      0.18 NLPi:   52 Dpth:    0 Lvs:    3 Obj: -2.81E+06 Gap:       NaN
Iter:     2 I:  0 Tm:      0.03 NLPi:   10 Dpth:    1 Lvs:    4 Obj: -2.44E+06 Gap:       NaN
Iter:     3 I:  0 Tm:      0.02 NLPi:   10 Dpth:    1 Lvs:    6 Obj: -2.66E+06 Gap:       NaN
Iter:     4 I:  0 Tm:     -0.00 NLPi:    7 Dpth:    1 Lvs:    8 Obj: -2.79E+06 Gap:       NaN
Iter:     5 I:  0 Tm:      0.03 NLPi:   15 Dpth:    2 Lvs:    9 Obj: -2.73E+06 Gap:       NaN
Iter:     6 I:  0 Tm:      0.04 NLPi:   19 Dpth:    2 Lvs:   11 Obj: -2.79E+06 Gap:       NaN
Iter:     7 I:  0 Tm:      0.02 NLPi:    3 Dpth:    2 Lvs:   12 Obj: -2.74E+06 Gap:       NaN
Iter:     8 I:  0 Tm:      0.03 NLPi:   16 Dpth:    3 Lvs:   14 Obj: -2.78E+06 Gap:       NaN
Iter:     9 I:  0 Tm:      0.00 NLPi:    3 Dpth:    3 Lvs:   16 Obj: -2.58E+06 Gap:       NaN
Iter:    10 I:  0 Tm:      0.00 NLPi:    3 Dpth:    3 Lvs:   18 Obj: -2.78E+06 Gap:       NaN
Iter:    11 I:  0 Tm:      0.04 NLPi:   16 Dpth:    4 Lvs:   19 Obj: -2.75E+06 Gap:       NaN
Iter:    12 I:  0 Tm:      0.00 NLPi:    5 Dpth:    4 Lvs:   20 Obj: -2.69E+06 Gap:       NaN
Iter:    13 I:  0 Tm:      0.02 NLPi:    5 Dpth:    4 Lvs:   21 Obj: -2.75E+06 Gap:       NaN
--Integer Solution:  -2.78E+06 Lowest Leaf:  -2.78E+06 Gap:   0.00E+00
Iter:    14 I:  0 Tm:      0.00 NLPi:    3 Dpth:    4 Lvs:   21 Obj: -2.78E+06 Gap:  0.00E+00
 Successful solution
 
 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :  0.45339999999999997 sec
 Objective      :  -2782980.2999999993
 Successful solution
 ---------------------------------------------------

Answer 2

在另一个答案的基础上，当将DCOST上的RegenRate1重新定义为MV时，有一个更好的解决方案(虽然计算时间为1.16秒，而不是0.45秒)。

​

​

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

m=GEKKO(remote=False)           #instantiate GEKKO object
m.time=np.linspace(0,23,24)     #define timespace

#Define production flow rates, when on, for each train
Train1_F=m.Var(value=300,lb=100,ub=350,name='Train1F')  #Train 1 Flow

#Helper variable representing regeneration rate
RegenRate1=m.MV(lb=0,ub=5800,name='RegenRate1')  #RegenRate1
RegenRate1.STATUS=1
RegenRate1.DCOST =0.01

#define regen state as binary variable --  1 is online, 0 is regeneration
Train1_On=m.MV(value=1,lb=0,ub=1,integer=True,name='Train1 Online')
Train1_On.STATUS=1

#define production variable as product of regen status and production rate - i.e 0 when on regen
Train1_Prod=m.Intermediate(Train1_F*Train1_On,name='Train1 Prod')

#maintain production accumulation per train,  variable per train,random starting points
Train1_Acc=m.Var(value=2000,lb=0,ub=5800,name='Train1_Acc')
RegenRate1_Off = m.Intermediate(RegenRate1*(1-Train1_On))
RegenRate1_On = m.Intermediate(RegenRate1*Train1_On)

#define accumulation as function of production - reset when regen is on
m.Equation(Train1_Acc.dt()==4*Train1_Prod - RegenRate1_Off)

#maximise production
Total_production=m.Intermediate(Train1_Prod)

#Attempts to Maximise production with as little as possible regen-cycles
m.Maximize(10*Total_production)
m.Maximize(100*RegenRate1_Off)
m.Minimize(RegenRate1_On)
m.Minimize(0.1*(1-Train1_On))

m.options.IMODE=6   #control mode
m.options.SOLVER=1  #select APOPT, suitable for Mixed Integer Programming

m.solver_options = ['minlp_gap_tol 1.0e-3',\
                    'minlp_maximum_iterations 10000',\
                    'minlp_max_iter_with_int_sol 500']

m.options.NODES=2
m.solve()

plt.figure(figsize=(12,5))
plt.subplot(2,1,1)
plt.plot(m.time,Train1_Acc,'b.-',label='Acc')
plt.plot(m.time,RegenRate1,'g.-',label='Regen')
plt.legend(); plt.grid()
plt.subplot(2,1,2)
plt.step(m.time,Train1_On,'r.-',label='Train On')
plt.legend(); plt.grid(); plt.xlabel('Time')
plt.show()