I have coded a parameter estimation of a 5-state Resistor-Capacitor (RC) model of a large building through Python GEKKO using 5 differential equations in MHE mode (IMODE 5). The code works fine but the problem is that the computation time becomes too high when the dataset is more than 1 week with 5 mins timesteps (34 parameters to estimate). More than 70% of the time is spent in "Dynamic Init. C" rather than in the solver itself (IPOPT) so I was wondering if someone knows what that "Dynamic Init. C" really means and how I could possibly speed up my model.
The model code is the following:
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
It is then applied to a for loop to verify the performance of this model for prediction, doing a new parameter estimation every 24 hours:
#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:]))
The prediction function is just the same model that uses the trained parameters to make predictions for the temperatures of the 5 building zones (IMODE 4). Almost all the time is spent computing the parameter estimation. The training computation time at every iteration is: [ 17.77, 34.25, 61.95, 97.31, 152.96, 215.53, 319.42, 463.88, 548.74, 691.37, 946.15, 1144.1 ] Last iteration takes around 19 minutes and it uses a 2 weeks dataset. The report of this last iteration, with the computation times included (DIAGLEVEL=1), is the following (i had to cut some parts because of characters limit):
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" takes 833 seconds on 1141. So that is definitely my bottleneck.
How I can improve it? I already tried to change IMODE or SOLVER but without success. I'll have to compute also an optimization (that will use also this model) for my MPC so it's quite important for me to reduce computation times.
Thanks to everyone that will answer to this post. If more data is needed to reproduce the code I could provide the rest of the training data (I would have to provide lot of code and some txt files)
Impressive application! Here is a summary of the model size for the nonlinear differential and algebraic equations once they are fully discretized:
Number of state variables: 179746
Number of total equations: - 179712
Number of slack variables: - 0
---------------------------------------
Degrees of freedom : 34
IPOPT requires 156.57 sec to solve while the model compilation requires 984.60 sec. The model building phases A, B, and C are the translation of the model into byte-code with automatic differentiation. The time grows with a longer time horizon because more instances of the model collocation equations are added with the longer horizon.
Without this byte-code compile, the model evaluation speed would be much slower if evaluated in Python. Here is a summary of how long the model evaluations took once they were compiled:
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)
There is a feature request in GitHub to store the compiled model so that it only needs to be built once. There is some functionality to keep the compiled model with the REPLAY option, but it isn't ready for Gekko MHE yet.
Here are a few recommendations to increase the speed with the current version of Gekko:
WMEAS (weight on measurements) and WMODEL (weight on prior model prediction). A higher WMODEL allows a shorter time horizon with similar results as you run the MHE over more data.m.options.TIME_SHIFT. Use 12 for a 1 hour time shift with 5 min sampled data.m.MV() to m.Param() because the additional MV equations are not needed for this application where it is only input data.Please post follow-up comments based on anything that helps or if you still need additional suggestions.