我使用 MHE 模式 (IMODE 5) 下的 5 个微分方程，通过 Python GEKKO 编写了大型建筑物 5 状态电阻电容 (RC) 模型的参数估计。 该代码工作正常，但问题是，当数据集超过 1 周、时间步长为 5 分钟(需要估计 34 个参数)时，计算时间会变得太长。超过 70% 的时间花在“Dynamic Init.C”而不是求解器本身(IPOPT)上，所以我想知道是否有人知道“Dynamic Init.C”的真正含义以及我如何才能加快我的速度型号。
模型代码如下:
def train (starting_datetime,finishing_datetime,i,i_max):
# Import or generate data
dates = pd.date_range(starting_datetime, finishing_datetime,freq=interp_timestep ,tz="US/Eastern")
df, location, solar_position, Q_solar_df = weather_data (starting_datetime, finishing_datetime)
U_train = Power_delivered[starting_datetime:finishing_datetime].values
X_train = Temp[starting_datetime:finishing_datetime].values
W_cols = ["T_ext", "Tground_WS", "Tground_NE", "Tground_G",
"Q_sol_WS", "Q_sol_NE", "Q_sol_WS_eta", "Q_sol_NE_eta",
"Q_int_WS", "Q_int_NE", "Q_int_G", "Q_int_WS_eta", "Q_int_NE_eta", "wind_sq"] # exogenous_columns
W_train = df[W_cols].values#.iloc[:-1].values
time_array = np.arange(0,len(dates)*timestep,timestep)
# Create GEKKO Model
m = GEKKO(remote=False)
m.time = time_array
# use previous solution as initial values after 5 iterations, before physics based values
if i > 5:
Uext_WS_initial = params[i-3]["Uext_WS"]
Uext_NE_initial = params[i-3]["Uext_NE"]
Uext_G_initial = params[i-3]["Uext_G"]
Uext_WS_eta_initial = params[i-3]["Uext_WS_eta"]
Uext_NE_eta_initial = params[i-3]["Uext_NE_eta"]
Ug_WS_initial = params[i-3]["Ug_WS"]
Ug_NE_initial = params[i-3]["Ug_NE"]
Ug_G_initial = params[i-3]["Ug_G"]
C_WS_initial = params[i-3]["C_WS"]
C_NE_initial = params[i-3]["C_NE"]
C_G_initial = params[i-3]["C_G"]
C_WS_eta_initial = params[i-3]["C_WS_eta"]
C_NE_eta_initial = params[i-3]["C_NE_eta"]
U_WS_initial = params[i-3]["U_WS"]
U_NE_initial = params[i-3]["U_NE"]
U_rdc_initial = params[i-3]["U_rdc"]
U_eta_initial = params[i-3]["U_eta"]
U_rdc_gym_initial = params[i-3]["U_rdc_gym"]
U_eta_gym_initial = params[i-3]["U_eta_gym"]
sigma_WS_initial = params[i-3]["sigma_WS"]
sigma_NE_initial = params[i-3]["sigma_NE"]
sigma_G_initial = params[i-3]["sigma_G"]
sigma_WS_eta_initial = params[i-3]["sigma_WS_eta"]
sigma_NE_eta_initial = params[i-3]["sigma_NE_eta"]
lambda_WS_initial = params[i-3]["lambda_WS"]
lambda_NE_initial = params[i-3]["lambda_NE"]
lambda_G_initial = params[i-3]["lambda_G"]
lambda_WS_eta_initial = params[i-3]["lambda_WS_eta"]
lambda_NE_eta_initial = params[i-3]["lambda_NE_eta"]
alpha_WS_initial = params[i-3]["alpha_WS"]
alpha_NE_initial = params[i-3]["alpha_NE"]
alpha_G_initial = params[i-3]["alpha_G"]
alpha_WS_eta_initial = params[i-3]["alpha_WS_eta"]
alpha_NE_eta_initial = params[i-3]["alpha_NE_eta"]
else:
Uext_WS_initial = Uext_WS_physics
Uext_NE_initial = Uext_NE_physics
Uext_G_initial = Uext_G_physics
Uext_WS_eta_initial = Uext_WS_eta_physics
Uext_NE_eta_initial = Uext_NE_eta_physics
Ug_WS_initial = Ug_WS_physics
Ug_NE_initial = Ug_NE_physics
Ug_G_initial = Ug_G_physics
C_WS_initial = C_WS_physics
C_NE_initial = C_NE_physics
C_G_initial = C_G_physics
C_WS_eta_initial = C_WS_eta_physics
C_NE_eta_initial = C_NE_eta_physics
U_WS_initial = U_WS_physics
U_NE_initial = U_NE_physics
U_rdc_initial = U_rdc_physics
U_eta_initial = U_eta_physics
U_rdc_gym_initial = U_rdc_gym_physics
U_eta_gym_initial = U_eta_gym_physics
sigma_WS_initial = 1
sigma_NE_initial = 1
sigma_G_initial = 1
sigma_WS_eta_initial = 1
sigma_NE_eta_initial = 1
lambda_WS_initial = 1
lambda_NE_initial = 1
lambda_G_initial = 1
lambda_WS_eta_initial = 1
lambda_NE_eta_initial = 1
alpha_WS_initial = alpha_WS_physics
alpha_NE_initial = alpha_NE_physics
alpha_G_initial = alpha_G_physics
alpha_WS_eta_initial = alpha_WS_eta_physics
alpha_NE_eta_initial = alpha_NE_eta_physics
# Parameters to Estimate
Uext_WS = m.FV(value=Uext_WS_initial,lb=0,ub=1.5) # kW / W
Uext_NE = m.FV(value=Uext_NE_initial,lb=0,ub=1.5) # kW / W
Uext_G = m.FV(value=Uext_G_initial,lb=0,ub=1.5) # kW / W
Uext_WS_eta = m.FV(value=Uext_WS_eta_initial,lb=0,ub=1.5) # kW / W
Uext_NE_eta = m.FV(value=Uext_NE_eta_initial,lb=0,ub=1.5) # kW / W
Ug_WS = m.FV(value=Ug_WS_initial,lb=0,ub=2) # kW / W
Ug_NE = m.FV(value=Ug_NE_initial,lb=0,ub=2) # kW / W
Ug_G = m.FV(value=Ug_G_initial,lb=0,ub=2) # kW / W
C_WS = m.FV(value=C_WS_initial,lb=0,ub=2) # GJ / K
C_NE = m.FV(value=C_NE_initial,lb=0,ub=2) # GJ / K
C_G = m.FV(value=C_G_initial,lb=0,ub=2) # GJ / K
C_WS_eta = m.FV(value=C_WS_eta_initial,lb=0,ub=2) # GJ / K
C_NE_eta = m.FV(value=C_NE_eta_initial,lb=0,ub=2) # GJ / K
U_WS = m.FV(value=U_WS_initial,lb=0,ub=2) # GJ / K
U_NE = m.FV(value=U_NE_initial,lb=0,ub=2) # GJ / K
U_rdc = m.FV(value=U_rdc_initial,lb=0,ub=1) # GJ / K
U_eta = m.FV(value=U_eta_initial,lb=0,ub=1) # GJ / K
U_rdc_gym = m.FV(value=U_rdc_gym_initial,lb=0,ub=1) # GJ / K
U_eta_gym = m.FV(value=U_eta_gym_initial,lb=0,ub=1) # GJ / K
sigma_WS = m.FV(value=sigma_WS_initial,lb=0,ub=10)
sigma_NE = m.FV(value=sigma_NE_initial,lb=0,ub=10)
sigma_G = m.FV(value=sigma_G_initial,lb=0,ub=10)
sigma_WS_eta = m.FV(value=sigma_WS_eta_initial,lb=0,ub=10)
sigma_NE_eta = m.FV(value=sigma_NE_eta_initial,lb=0,ub=10)
lambda_WS = m.FV(value=lambda_WS_initial,lb=0,ub=10)
lambda_NE = m.FV(value=lambda_NE_initial,lb=0,ub=10)
lambda_G = m.FV(value=lambda_G_initial,lb=0,ub=10)
lambda_WS_eta = m.FV(value=lambda_WS_eta_initial,lb=0,ub=10)
lambda_NE_eta = m.FV(value=lambda_NE_eta_initial,lb=0,ub=10)
alpha_WS = m.FV(value=alpha_WS_initial,lb=0,ub=0.3)
alpha_NE = m.FV(value=alpha_NE_initial,lb=0,ub=0.3)
alpha_G = m.FV(value=alpha_G_initial,lb=0,ub=0.3)
alpha_WS_eta = m.FV(value=alpha_WS_eta_initial,lb=0,ub=0.3)
alpha_NE_eta = m.FV(value=alpha_NE_eta_initial,lb=0,ub=0.3)
# STATUS=1 allows solver to adjust parameter
Uext_WS.STATUS = 1
Uext_NE.STATUS = 1
Uext_G.STATUS = 1
Uext_WS_eta.STATUS = 1
Uext_NE_eta.STATUS = 1
Ug_WS.STATUS = 1
Ug_NE.STATUS = 1
Ug_G.STATUS = 1
C_WS.STATUS = 1
C_NE.STATUS = 1
C_G.STATUS = 1
C_WS_eta.STATUS = 1
C_NE_eta.STATUS = 1
U_WS.STATUS = 1
U_NE.STATUS = 1
U_rdc.STATUS = 1
U_eta.STATUS = 1
U_rdc_gym.STATUS = 1
U_eta_gym.STATUS = 1
sigma_WS.STATUS = 1
sigma_NE.STATUS = 1
sigma_G.STATUS = 1
sigma_WS_eta.STATUS = 1
sigma_NE_eta.STATUS = 1
lambda_WS.STATUS = 1
lambda_NE.STATUS = 1
lambda_G.STATUS = 1
lambda_WS_eta.STATUS = 1
lambda_NE_eta.STATUS = 1
alpha_WS.STATUS = 1
alpha_NE.STATUS = 1
alpha_G.STATUS = 1
alpha_WS_eta.STATUS = 1
alpha_NE_eta.STATUS = 1
# Measured inputs
Q_WS = m.MV(value=U_train[:,0])
Q_NE = m.MV(value=U_train[:,1])
Q_G = m.MV(value=U_train[:,2])
Q_WS_eta = m.MV(value=U_train[:,3])
Q_NE_eta = m.MV(value=U_train[:,4])
T_ext = m.MV(value=W_train[:,0])
Tground_WS = m.MV(value=W_train[:,1])
Tground_NE = m.MV(value=W_train[:,2])
Tground_G = m.MV(value=W_train[:,3])
Q_sol_WS = m.MV(value=W_train[:,4])
Q_sol_NE = m.MV(value=W_train[:,5])
Q_sol_WS_eta = m.MV(value=W_train[:,6])
Q_sol_NE_eta = m.MV(value=W_train[:,7])
Q_int_WS = m.MV(value=W_train[:,8])
Q_int_NE = m.MV(value=W_train[:,9])
Q_int_G = m.MV(value=W_train[:,10])
Q_int_WS_eta = m.MV(value=W_train[:,11])
Q_int_NE_eta = m.MV(value=W_train[:,12])
wind_sq = m.MV(value=W_train[:,13])
# State variables
T_WS = m.CV(value=X_train[:,0])
T_NE = m.CV(value=X_train[:,1])
T_G = m.CV(value=X_train[:,2])
T_WS_eta = m.CV(value=X_train[:,3])
T_NE_eta = m.CV(value=X_train[:,4])
T_WS.FSTATUS = 1 # minimize fstatus * (meas-pred)^2
T_NE.FSTATUS = 1 # minimize fstatus * (meas-pred)^2
T_G.FSTATUS = 1 # minimize fstatus * (meas-pred)^2
T_WS_eta.FSTATUS = 1 # minimize fstatus * (meas-pred)^2
T_NE_eta.FSTATUS = 1 # minimize fstatus * (meas-pred)^2
# Semi-fundamental correlations (energy balances)
m.Equation(C_WS *1e6* T_WS.dt() == (Q_WS + lambda_WS * Q_int_WS + sigma_WS * Q_sol_WS - alpha_WS * wind_sq \
+ Uext_WS*(T_ext - T_WS) + Ug_WS*(Tground_WS - T_WS) \
+ U_rdc*(T_NE - T_WS) + U_rdc_gym*(T_G - T_WS) + U_WS*(T_WS_eta - T_WS) ))
m.Equation(C_NE *1e6* T_NE.dt() == (Q_NE + lambda_NE * Q_int_NE + sigma_NE * Q_sol_NE - alpha_NE * wind_sq \
+ Uext_NE*(T_ext - T_NE) + Ug_NE*(Tground_NE - T_NE) \
+ U_rdc*(T_WS - T_NE) + U_NE*(T_NE_eta - T_NE) ))
m.Equation(C_G *1e6* T_G.dt() == (Q_G + lambda_G * Q_int_G + sigma_G * Q_sol_WS - alpha_G * wind_sq \
+ Uext_G*(T_ext - T_G) + Ug_G*(Tground_G - T_G) \
+ U_rdc_gym*(T_WS - T_G) + U_eta_gym*(T_WS_eta - T_G) ))
m.Equation(C_WS_eta *1e6* T_WS_eta.dt() == (Q_WS_eta + lambda_WS_eta * Q_int_WS_eta + sigma_WS_eta * Q_sol_WS_eta - alpha_WS_eta * wind_sq \
+ Uext_WS_eta*(T_ext - T_WS_eta) \
+ U_eta*(T_NE_eta - T_WS_eta) + U_eta_gym*(T_G - T_WS_eta) + U_WS*(T_WS - T_WS_eta) ))
m.Equation(C_NE_eta *1e6* T_NE_eta.dt() == (Q_NE_eta + lambda_NE_eta * Q_int_NE_eta + sigma_NE_eta * Q_sol_NE_eta - alpha_NE_eta * wind_sq \
+ Uext_NE_eta*(T_ext - T_NE_eta) \
+ U_eta*(T_WS_eta - T_NE_eta) + U_NE*(T_NE - T_NE_eta) ))
# Options
m.options.IMODE = 5 # MHE
m.options.EV_TYPE = 2 # Objective type
m.options.NODES = 1 # Collocation nodes
m.options.SOLVER = 3 # IPOPT
# Solve
if (i == i_max and i != 1):
m.options.DIAGLEVEL = 1
m.solve(disp=True)
else:
m.solve(disp=False)
train_results = {"Uext_WS": Uext_WS.value[0], "Uext_NE": Uext_NE.value[0], "Uext_G": Uext_G.value[0], "Uext_WS_eta": Uext_WS_eta.value[0],
"Uext_NE_eta": Uext_NE_eta.value[0], "Ug_WS": Ug_WS.value[0], "Ug_NE": Ug_NE.value[0], "Ug_G": Ug_G.value[0],
"C_WS": C_WS.value[0], "C_NE": C_NE.value[0], "C_G": C_G.value[0], "C_WS_eta": C_WS_eta.value[0],
"C_NE_eta": C_NE_eta.value[0], "U_WS": U_WS.value[0], "U_NE": U_NE.value[0], "U_rdc": U_rdc.value[0],
"U_eta": U_eta.value[0], "U_rdc_gym": U_rdc_gym.value[0], "U_eta_gym": U_eta_gym.value[0],
"sigma_WS": sigma_WS.value[0], "sigma_NE": sigma_NE.value[0], "sigma_G": sigma_G.value[0],
"sigma_WS_eta": sigma_WS_eta.value[0], "sigma_NE_eta": sigma_NE_eta.value[0], "lambda_WS": lambda_WS.value[0],
"lambda_NE": lambda_NE.value[0], "lambda_G": lambda_G.value[0], "lambda_WS_eta": lambda_WS_eta.value[0],
"lambda_NE_eta": lambda_NE_eta.value[0], "alpha_WS": alpha_WS.value[0], "alpha_NE": alpha_NE.value[0],
"alpha_G": alpha_G.value[0], "alpha_WS_eta": alpha_WS_eta.value[0], "alpha_NE_eta": alpha_NE_eta.value[0]}
return train_results
然后将其应用于 for 循环以验证该模型的预测性能，每 24 小时进行一次新的参数估计:
#first iteration and initialization
starting_datetime = pd.Timestamp("2015-02-01 00:00", tz="US/Eastern")
n_days_pred = 1
i = 1
i_max = 1
n_timesteps = int(n_days_pred*24*3600 / timestep)
finishing_datetime = starting_datetime + timedelta(hours = (24*n_days_pred))
train_results = train (starting_datetime,finishing_datetime,i,i_max)
finishing_datetime = starting_datetime + timedelta(hours = (24*n_days_pred*2))
X_pred = prediction (train_results, starting_datetime,finishing_datetime)
#the rest of the iterations
training_time = []
params=[]
i_max = 13
for i in range (2,i_max+1):
finishing_datetime = starting_datetime + timedelta(hours = (24*n_days_pred*i))
start_train = time.time()
train_results = train (starting_datetime,finishing_datetime,i,i_max)
print(train_results)
end_train = time.time()
training_time.append(end_train - start_train)
print(i)
print(training_time)
params.append(train_results)
starting_datetime_pred = finishing_datetime
finishing_datetime_pred = finishing_datetime + timedelta(hours = (24*n_days_pred*1))
X_pred_new = prediction (train_results, starting_datetime_pred,finishing_datetime_pred)
X_pred = np.concatenate((X_pred,X_pred_new[-n_timesteps:]))
预测函数与使用训练参数来预测 5 个建筑区域的温度的模型相同 (IMODE 4)。 几乎所有的时间都花在计算参数估计上。 每次迭代的训练计算时间为: [ 17.77、34.25、61.95、97.31、152.96、215.53、319.42、 463.88、548.74、691.37、946.15、1144.1] 最后一次迭代大约需要 19 分钟，并且使用 2 周的数据集。最后一次迭代的报告，包括计算时间(DIAGLEVEL=1)，如下(由于字符限制，我不得不删除一些部分):
Called files( 35 )
READ info FILE FOR VARIABLE DEFINITION: gk_model116.info
initialization step: 0
Parsing model file gk_model116.apm
Read model file (sec): 0.013000000000000012
Initialize constants (sec): 0.
Determine model size (sec): 0.006000000000000005
Allocate memory (sec): 0.
Parse and store model (sec): 0.0030000000000000027
--------- APM Model Size ------------
Each time step contains
Objects : 0
Constants : 0
Variables : 58
Intermediates: 0
Connections : 0
Equations : 5
Residuals : 5
Error checking (sec): 0.0010000000000000009
Compile equations (sec): 0.002999999999999975
Check for uninitialized intermediates (sec): 0.
------------------------------------------------------
Total Parse Time (sec): 0.027000000000000007
initialization step: 1
initialization step: 2
initialization step: 3
Called files( 8 )
Files(8): File Read ss.t0 F
files: ss.t0 does not exist
Called files( 31 )
READ info FILE FOR PROBLEM DEFINITION: gk_model116.info
initialization step: 4
Called files( 31 )
READ info FILE FOR PROBLEM DEFINITION: gk_model116.info
WARNING: restart file not found - est.t0
Called files( 51 )
Read DBS File defaults.dbs
files: defaults.dbs does not exist
Called files( 51 )
Read DBS File gk_model116.dbs
files: gk_model116.dbs does not exist
Called files( 51 )
Read DBS File measurements.dbs
Called files( 51 )
Read DBS File overrides.dbs
files: overrides.dbs does not exist
Number of state variables: 179746
Number of total equations: - 179712
Number of slack variables: - 0
---------------------------------------
Degrees of freedom : 34
This is Ipopt version 3.10.2, running with linear solver mumps.
Number of nonzeros in equality constraint Jacobian...: 501691
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 183456
Total number of variables............................: 179746
variables with only lower bounds: 0
variables with lower and upper bounds: 34
variables with only upper bounds: 0
Total number of equality constraints.................: 179712
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
Number of Iterations....: 48
(scaled) (unscaled)
Objective...............: 1.4797975081572714e+005 1.4797975081572714e+005
Dual infeasibility......: 1.4113310324500653e-009 1.4113310324500653e-009
Constraint violation....: 1.6653345369377348e-015 4.2632564145606011e-014
Complementarity.........: 1.0942671862849926e-008 1.0942671862849926e-008
Overall NLP error.......: 7.3565689235462651e-010 1.0942671862849926e-008
Number of objective function evaluations = 56
Number of objective gradient evaluations = 49
Number of equality constraint evaluations = 56
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 49
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 48
Total CPU secs in IPOPT (w/o function evaluations) = 65.561
Total CPU secs in NLP function evaluations = 80.763
EXIT: Optimal Solution Found.
The solution was found.
The final value of the objective function is 147979.75081572714
---------------------------------------------------
Solver : IPOPT (v3.12)
Solution time : 156.5746999999999 sec
Objective : 147979.75081572714
Successful solution
---------------------------------------------------
Called files( 2 )
Called files( 65 )
Called files( 66 )
Called files( 83 )
Called files( 81 )
Called files( 52 )
WRITE dbs FILE
Called files( 56 )
WRITE json FILE
Timer # 1 1141.17/ 1 = 1141.17 Total system time
Timer # 2 156.57/ 1 = 156.57 Total solve time
Timer # 3 11.16/ 56 = 0.20 Objective Calc: apm_p
Timer # 4 11.36/ 50 = 0.23 Objective Grad: apm_g
Timer # 5 9.31/ 56 = 0.17 Constraint Calc: apm_c
Timer # 6 0.00/ 1 = 0.00 Sparsity: apm_s
Timer # 7 7.96/ 50 = 0.16 1st Deriv #1: apm_a1
Timer # 8 0.00/ 0 = 0.00 1st Deriv #2: apm_a2
Timer # 9 0.91/ 7488 = 0.00 Custom Init: apm_custom_init
Timer # 10 0.19/ 7488 = 0.00 Mode: apm_node_res::case 0
Timer # 11 0.12/ 22464 = 0.00 Mode: apm_node_res::case 1
Timer # 12 0.46/ 7488 = 0.00 Mode: apm_node_res::case 2
Timer # 13 0.06/ 14976 = 0.00 Mode: apm_node_res::case 3
Timer # 14 4.80/ 883584 = 0.00 Mode: apm_node_res::case 4
Timer # 15 11.88/ 748800 = 0.00 Mode: apm_node_res::case 5
Timer # 16 10.17/ 366912 = 0.00 Mode: apm_node_res::case 6
Timer # 17 2.34/ 50 = 0.05 Base 1st Deriv: apm_jacobian
Timer # 18 0.00/ 0 = 0.00 Base 1st Deriv: apm_condensed_jacobian
Timer # 19 0.06/ 1 = 0.06 Non-zeros: apm_nnz
Timer # 20 0.00/ 0 = 0.00 Count: Division by zero
Timer # 21 0.00/ 0 = 0.00 Count: Argument of LOG10 negative
Timer # 22 0.00/ 0 = 0.00 Count: Argument of LOG negative
Timer # 23 0.00/ 0 = 0.00 Count: Argument of SQRT negative
Timer # 24 0.00/ 0 = 0.00 Count: Argument of ASIN illegal
Timer # 25 0.00/ 0 = 0.00 Count: Argument of ACOS illegal
Timer # 26 0.05/ 1 = 0.05 Extract sparsity: apm_sparsity
Timer # 27 0.41/ 17 = 0.02 Variable ordering: apm_var_order
Timer # 28 0.00/ 0 = 0.00 Condensed sparsity
Timer # 29 9.32/ 1 = 9.32 Hessian Non-zeros
Timer # 30 0.54/ 4 = 0.14 Differentials
Timer # 31 7.63/ 48 = 0.16 Hessian Calculation
Timer # 32 1.89/ 49 = 0.04 Extract Hessian
Timer # 33 0.19/ 1 = 0.19 Base 1st Deriv: apm_jac_order
Timer # 34 0.24/ 1 = 0.24 Solver Setup
Timer # 35 66.11/ 1 = 66.11 Solver Solution
Timer # 36 2.30/ 68 = 0.03 Number of Variables
Timer # 37 0.14/ 5 = 0.03 Number of Equations
Timer # 38 1.42/ 19 = 0.07 File Read/Write
Timer # 39 0.01/ 1 = 0.01 Dynamic Init A
Timer # 40 138.16/ 1 = 138.16 Dynamic Init B
Timer # 41 833.15/ 1 = 833.15 Dynamic Init C
Timer # 42 0.01/ 1 = 0.01 Init: Read APM File
Timer # 43 0.00/ 1 = 0.00 Init: Parse Constants
Timer # 44 0.01/ 1 = 0.01 Init: Model Sizing
Timer # 45 0.00/ 1 = 0.00 Init: Allocate Memory
Timer # 46 0.00/ 1 = 0.00 Init: Parse Model
Timer # 47 0.00/ 1 = 0.00 Init: Check for Duplicates
Timer # 48 0.00/ 1 = 0.00 Init: Compile Equations
Timer # 49 0.00/ 1 = 0.00 Init: Check Uninitialized
Timer # 50 0.00/ 7590 = 0.00 Evaluate Expression Once
Timer # 51 0.00/ 0 = 0.00 Sensitivity Analysis: LU Factorization
Timer # 52 0.00/ 0 = 0.00 Sensitivity Analysis: Gauss Elimination
Timer # 53 0.00/ 0 = 0.00 Sensitivity Analysis: Total Time
“Dynamic Init.C”在 1141 上花费了 833 秒。所以这绝对是我的瓶颈。
我该如何改进它？ 我已经尝试更改 IMODE 或 SOLVER 但没有成功。 我还必须为我的 MPC 计算优化(也将使用此模型)，因此减少计算时间对我来说非常重要。
感谢所有回答这篇文章的人。如果需要更多数据来重现代码，我可以提供其余的训练数据(我必须提供大量代码和一些 txt 文件)
最佳答案
令人印象深刻的应用!以下是非线性微分方程和代数方程完全离散化后的模型大小的摘要:
Number of state variables: 179746
Number of total equations: - 179712
Number of slack variables: - 0
---------------------------------------
Degrees of freedom : 34
IPOPT 需要 156.57 秒来求解，而模型编译需要 984.60 秒。模型构建阶段 A、B 和 C 是将模型通过自动微分转换为字节码。时间随着时间范围的延长而增长，因为随着时间范围的延长，会添加更多的模型配置方程实例。
如果没有这种字节码编译，如果用 Python 进行评估，模型评估速度会慢很多。以下是模型评估编译后所需时间的摘要:
Timer # 14 4.80 sec for 883584 times (eqn Residuals)
Timer # 15 11.88 sec for 748800 times (sparse Jacobian)
Timer # 16 10.17 sec for 366912 times (sparse Hessian)
有一个feature request in GitHub存储编译后的模型，使其只需要构建一次。有一些功能可以将编译后的模型保留为 REPLAY option ，但尚未为 Gekko MHE 做好准备。
以下是一些提高当前版本 Gekko 速度的建议:
- 在一周内解决问题时，将时间范围的稀疏性增加到 10-60 分钟时间步
- 包括一个遗忘因子来近似过去的数据和更长的时间范围的影响。平衡参数是 WMEAS(测量权重)和 WMODEL(先前模型预测权重)。当您对更多数据运行 MHE 时，较高的
WMODEL可以缩短时间范围并获得相似的结果。 - 每小时计算一次新参数，并使用
m.options.TIME_SHIFT控制数据和模型初始条件的时移。使用12进行 1 小时的时移和 5 分钟的采样数据。 - 将
m.MV()更改为m.Param()，因为此应用程序不需要额外的MV方程仅输入数据。
请根据任何有帮助的内容或如果您仍需要其他建议来发表后续评论。
关于Python Gekko 参数估计 : too long computational time because of "Dynamic Init C"，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/76807167/