英文:
Grey-box parameter estimation with GEKKO (@error: Insufficient Model)
问题
我使用GEKKO来估算热储罐的参数。我将储罐分成了20个等温层，并为每个等温层编写了能量平衡方程。未知参数包括相邻层之间的热导率和与环境的热损失系数。以下是您的代码的翻译：
# 导入必要的库
import numpy as np
import pandas as pd
from gekko import GEKKO
# 从CSV文件中读取数据
df = pd.read_csv('./dataset/charge_2.csv')
# 创建GEKKO模型
m = GEKKO()
# 参数
rho = m.Const(1000) # 水的密度 [kg/m^3]
r = m.Const(0.6) # 储罐的内半径
Ac = m.Const(np.pi * (r.value * r.value)) # 储罐的横截面积 [m^2]
Cp = m.Const(4186) # 水的比热容 [J/(kg.K)]
H = m.Const(1.79) # 储罐的高度 [m]
w = m.Const(H.value / 20) # 每个层的宽度 [m]
M = m.Const(rho.value * Ac.value * w.value) # 每个层的质量 [kg]
As = m.Const(2 * np.pi * r.value * w.value) # 每个层的表面积
Dt = m.Const(60) # 采样时间 [sec]
N = len(df.mass_flowrate_charge.values) # 时间步数
# 输入
m_dot = m.Array(m.Param, N) # 质量流率 [kg/s]
T_in = m.Array(m.Param, N) # 入口温度 [deg C]
Tamb_top = m.Array(m.Param, N) # 顶部环境温度 [deg C]
Tamb_mid = m.Array(m.Param, N) # 中部环境温度 [deg C]
Tamb_bot = m.Array(m.Param, N) # 底部环境温度 [deg C]
for i in range(N):
m_dot[i].value = df.mass_flowrate_charge.values[i]
T_in[i].value = df.Tin_charge.values[i]
Tamb_top[i].value = df.Tamb_top.values[i]
Tamb_bot[i].value = df.Tamb_bot.values[i]
Tamb_mid[i].value = df.Tamb_mid.values[i]
# 观测值
T_out = m.Array(m.Param, N) # 出口温度 [deg C]
T_top = m.Array(m.Param, N) # 顶部传感器温度 [deg C]
T_mid = m.Array(m.Param, N) # 中部传感器温度 [deg C]
T_bot = m.Array(m.Param, N) # 底部传感器温度 [deg C]
for i in range(N):
T_out[i].value = df.Tout_charge.values[i]
T_top[i].value = df.Ttop_m.values[i]
T_mid[i].value = df.Tmid_m.values[i]
T_bot[i].value = df.Tbot_m.values[i]
# 初始温度值对所有层
T0 = T_top[0].value
# 初始化用于存储结果的数组（层的温度）
T_L1, T_L2, T_L3, T_L4, T_L5, T_L6, T_L7, T_L8, T_L9, T_L10,\
T_L11, T_L12, T_L13, T_L14, T_L15, T_L16, T_L17, T_L18, T_L19, T_L20\
= m.Array(m.Var, N+1), m.Array(m.Var, N+1), m.Array(m.Var, N+1), \
m.Array(m.Var, N+1), m.Array(m.Var, N+1), m.Array(m.Var, N+1),\
m.Array(m.Var, N+1), m.Array(m.Var, N+1), m.Array(m.Var, N+1),\
m.Array(m.Var, N+1), m.Array(m.Var, N+1), m.Array(m.Var, N+1),\
m.Array(m.Var, N+1), m.Array(m.Var, N+1), m.Array(m.Var, N+1),\
m.Array(m.Var, N+1), m.Array(m.Var, N+1), m.Array(m.Var, N+1),\
m.Array(m.Var, N+1), m.Array(m.Var, N+1)
# 设置初始条件
T_L1[0].value, T_L2[0].value, T_L3[0].value, T_L4[0].value, T_L5[0].value,\
T_L6[0].value, T_L7[0].value, T_L8[0].value, T_L9[0].value, T_L10[0].value,\
T_L11[0].value, T_L12[0].value, T_L13[0].value, T_L14[0].value, T_L15[0].value,\
T_L16[0].value, T_L17[0].value, T_L18[0].value, T_L19[0].value, T_L20[0].value\
= T0, T0, T0, T0, T0, T0, T0, T0, T0, T0, T0, T0, T0, T0, T0, T0, T0, T0, T0, T0
# 可识别的参数
# k_{i,i+1} = 第i层与第i+1层之间的热导率 [W/(m.K)]。ka - 储罐的热系数 [W/m2]
k_12, k_23, k_34, k_45, k_56, k_67, k_78, k_89, k_910, k_1011, k_1112, k_1213,\
k_1314, k_1415, k_1516, k_1617, k_1718, k_1819, k_1920, ka = \
m.Var(0.6, lb=0, ub=100), m.Var(0.6, lb=0, ub=100), m.Var(0.6, lb=0, ub=100),\
m.Var(0.6, lb=
<details>
<summary>英文:</summary>
I am using GEKKO to estimate the parameters of a thermal storage tank. I divided the tank into 20 isothermal layers and have written the energy balance for each. The unknown parameters are the thermal conductivity between adjancent layers and the thermal loss coefficient to the environment. My code is as follows:
import numpy as np
import pandas as pd
from gekko import GEKKO
df = pd.read_csv('./dataset/charge_2.csv')
create GEKKO model
m = GEKKO()
Parameters
rho = m.Const(1000) # density of water [kg/m^3]
r = m.Const(0.6) # inner radius of the tank
Ac = m.Const(np.pi*(r.valuer.value)) # cross-sectional area of tank [m^2]
Cp = m.Const(4186) # specific heat capacity of water [J/(kg.K)]
H = m.Const(1.79) # height of the tank [m]
w = m.Const(H.value/20) # Width of each layer [m]
M = m.Const(rho.valueAc.valuew.value) # mass of each layer [kg]
As = m.Const(2np.pir.valuew.value) # surface area of each layer
Dt = m.Const(60) # Sample time [sec]
N = len(df.mass_flowrate_charge.values) # Number of time steps
Inputs
m_dot = m.Array(m.Param, N) # mass flow rate [kg/s]
T_in = m.Array(m.Param, N) # inlet temperature [deg C]
Tamb_top = m.Array(m.Param, N) # ambient temperature at top [deg C]
Tamb_mid = m.Array(m.Param, N) # ambient temperature at mid [deg C]
Tamb_bot = m.Array(m.Param, N) # ambient temperature at bottom [deg C]
for i in range(N):
m_dot[i].value = df.mass_flowrate_charge.values[i]
T_in[i].value = df.Tin_charge.values[i]
Tamb_top[i].value = df.Tamb_top.values[i]
Tamb_bot[i].value = df.Tamb_bot.values[i]
Tamb_mid[i].value = df.Tamb_mid.values[i]
Observations
T_out = m.Array(m.Param, N) # outlet temperature [deg C]
T_top = m.Array(m.Param, N) # top sensor temperature [deg C]
T_mid = m.Array(m.Param, N) # mid sensor temperature [deg C]
T_bot = m.Array(m.Param, N) # bottom sensor temperature [deg C]
for i in range(N):
T_out[i].value = df.Tout_charge.values[i]
T_top[i].value = df.Ttop_m.values[i]
T_mid[i].value = df.Tmid_m.values[i]
T_bot[i].value = df.Tbot_m.values[i]
Initial temperature for all layers
T0 = T_top[0].value
Initialize arrays to store results (temperature of layers)
T_L1, T_L2, T_L3, T_L4, T_L5, T_L6, T_L7, T_L8, T_L9, T_L10,
T_L11, T_L12, T_L13, T_L14, T_L15, T_L16, T_L17, T_L18, T_L19, T_L20
= m.Array(m.Var, N+1), m.Array(m.Var, N+1), m.Array(m.Var, N+1),
m.Array(m.Var, N+1), m.Array(m.Var, N+1), m.Array(m.Var, N+1),
m.Array(m.Var, N+1), m.Array(m.Var, N+1), m.Array(m.Var, N+1),
m.Array(m.Var, N+1), m.Array(m.Var, N+1), m.Array(m.Var, N+1),
m.Array(m.Var, N+1), m.Array(m.Var, N+1), m.Array(m.Var, N+1),
m.Array(m.Var, N+1), m.Array(m.Var, N+1), m.Array(m.Var, N+1),
m.Array(m.Var, N+1), m.Array(m.Var, N+1)
Set initial conditions
T_L1[0].value, T_L2[0].value, T_L3[0].value, T_L4[0].value, T_L5[0].value,
T_L6[0].value, T_L7[0].value, T_L8[0].value, T_L9[0].value, T_L10[0].value,
T_L11[0].value, T_L12[0].value, T_L13[0].value, T_L14[0].value, T_L15[0].value,
T_L16[0].value, T_L17[0].value, T_L18[0].value, T_L19[0].value, T_L20[0].value
= T0, T0, T0, T0, T0, T0, T0, T0, T0, T0, T0, T0, T0, T0, T0, T0, T0, T0, T0, T0
Identifiable parameters
k_{i,i+1} = thermal conductivity between layer i and i+1 [W/(m.K)]. ka - thermal
coefficient of tank [W/m2]
k_12, k_23, k_34, k_45, k_56, k_67, k_78, k_89, k_910, k_1011, k_1112, k_1213,
k_1314, k_1415, k_1516, k_1617, k_1718, k_1819, k_1920, ka =
m.Var(0.6, lb = 0, ub = 100), m.Var(0.6, lb = 0, ub = 100), m.Var(0.6, lb = 0, ub = 100),
m.Var(0.6, lb = 0, ub = 100), m.Var(0.6, lb = 0, ub = 100), m.Var(0.6, lb = 0, ub = 100),
m.Var(0.6, lb = 0, ub = 100), m.Var(0.6, lb = 0, ub = 100), m.Var(0.6, lb = 0, ub = 100),
m.Var(0.6, lb = 0, ub = 100), m.Var(0.6, lb = 0, ub = 100), m.Var(0.6, lb = 0, ub = 100),
m.Var(0.6, lb = 0, ub = 100), m.Var(0.6, lb = 0, ub = 100), m.Var(0.6, lb = 0, ub = 100),
m.Var(0.6, lb = 0, ub = 100), m.Var(0.6, lb = 0, ub = 100), m.Var(0.6, lb = 0, ub = 100),
m.Var(0.6, lb = 0, ub = 100), m.Var(5, lb = 0, ub = 100)\
Solve system of equations using Euler's method
for n in range(N):
# Energy Balance for each layer (conduction between layer + thermal losses
# to environment + heat transfer due to fluid flow in layers)
dT_L1_dt = m.Intermediate(((k_12*Ac)/w*(T_L2[n] - T_L1[n])
+ ka*(Tamb_top[n] - T_L1[n])*(As+Ac)
)/(M*Cp))
dT_L2_dt = m.Intermediate(((k_12*Ac)/w*(T_L1[n] - T_L2[n]) + (k_23*Ac)/w*(T_L3[n] - T_L2[n])
+ ka*(Tamb_top[n] - T_L2[n])*As
)/(M*Cp))
dT_L3_dt = m.Intermediate(((k_23*Ac)/w*(T_L2[n] - T_L3[n]) + (k_34*Ac)/w*(T_L4[n] - T_L3[n])
+ m_dot[n]*Cp*(T_in[n] - T_L3[n])
+ ka*(Tamb_top[n] - T_L3[n])*As
)/(M*Cp))
dT_L4_dt = m.Intermediate(((k_34*Ac)/w*(T_L3[n] - T_L4[n]) + (k_45*Ac)/w*(T_L5[n] - T_L4[n])
+ m_dot[n]*Cp*(T_L3[n] - T_L4[n])
+ ka*(Tamb_top[n] - T_L4[n])*As
)/(M*Cp))
dT_L5_dt = m.Intermediate(((k_45*Ac)/w*(T_L4[n] - T_L5[n]) + (k_56*Ac)/w*(T_L6[n] - T_L5[n])
+ m_dot[n]*Cp*(T_L4[n] - T_L5[n])
+ ka*(Tamb_top[n] - T_L5[n])
)/(M*Cp))
dT_L6_dt = m.Intermediate(((k_56*Ac)/w*(T_L5[n] - T_L6[n]) + (k_67*Ac)/w*(T_L7[n] - T_L6[n])
+ m_dot[n]*Cp*(T_L5[n] - T_L6[n])
+ ka*(Tamb_top[n] - T_L6[n])*As
)/(M*Cp))
dT_L7_dt = m.Intermediate(((k_67*Ac)/w*(T_L6[n] - T_L7[n]) + (k_78*Ac)/w*(T_L8[n] - T_L7[n])
+ m_dot[n]*Cp*(T_L6[n] - T_L7[n])
+ ka*(Tamb_top[n] - T_L7[n])*As
)/(M*Cp))
dT_L8_dt = m.Intermediate(((k_78*Ac)/w*(T_L7[n] - T_L8[n]) + (k_89*Ac)/w*(T_L9[n] - T_L8[n])
+ m_dot[n]*Cp*(T_L7[n] - T_L8[n])
+ ka*(Tamb_mid[n] - T_L8[n])*As
)/(M*Cp))
dT_L9_dt = m.Intermediate(((k_89*Ac)/w*(T_L8[n] - T_L9[n]) + (k_910*Ac)/w*(T_L10[n] - T_L9[n])
+ m_dot[n]*Cp*(T_L8[n] - T_L9[n])
+ ka*(Tamb_mid[n] - T_L9[n])*As
)/(M*Cp))
dT_L10_dt = m.Intermediate(((k_910*Ac)/w*(T_L9[n] - T_L10[n]) + (k_1011*Ac)/w*(T_L11[n] - T_L10[n])
+ m_dot[n]*Cp*(T_L9[n] - T_L10[n])
+ ka*(Tamb_mid[n] - T_L10[n])*As
)/(M*Cp))
dT_L11_dt = m.Intermediate(((k_1011*Ac)/w*(T_L10[n] - T_L11[n]) + (k_1112*Ac)/w*(T_L12[n] - T_L11[n])
+ m_dot[n]*Cp*(T_L10[n] - T_L11[n])
+ ka*(Tamb_mid[n] - T_L11[n])*As
)/(M*Cp))
dT_L12_dt = m.Intermediate(((k_1112*Ac)/w*(T_L11[n] - T_L12[n]) + (k_1213*Ac)/w*(T_L13[n] - T_L12[n])
+ m_dot[n]*Cp*(T_L11[n] - T_L12[n])
+ ka*(Tamb_mid[n] - T_L12[n])*As
)/(M*Cp))
dT_L13_dt = m.Intermediate(((k_1213*Ac)/w*(T_L12[n] - T_L13[n]) + (k_1314*Ac)/w*(T_L14[n] - T_L13[n])
+ m_dot[n]*Cp*(T_L12[n] - T_L13[n])
+ ka*(Tamb_mid[n] - T_L13[n])*As
)/(M*Cp))
dT_L14_dt = m.Intermediate(((k_1314*Ac)/w*(T_L13[n] - T_L14[n]) + (k_1415*Ac)/w*(T_L15[n] - T_L14[n])
+ m_dot[n]*Cp*(T_L13[n] - T_L14[n])
+ ka*(Tamb_bot[n] - T_L14[n])*As
)/(M*Cp))
dT_L15_dt = m.Intermediate(((k_1415*Ac)/w*(T_L14[n] - T_L15[n]) + (k_1516*Ac)/w*(T_L16[n] - T_L15[n])
+ m_dot[n]*Cp*(T_L14[n] - T_L15[n])
+ ka*(Tamb_bot[n] - T_L15[n])*As
)/(M*Cp))
dT_L16_dt = m.Intermediate(((k_1516*Ac)/w*(T_L15[n] - T_L16[n]) + (k_1617*Ac)/w*(T_L17[n] - T_L16[n])
+ m_dot[n]*Cp*(T_L15[n] - T_L16[n])
+ ka*(Tamb_bot[n] - T_L16[n])*As
)/(M*Cp))
dT_L17_dt = m.Intermediate(((k_1617*Ac)/w*(T_L16[n] - T_L17[n]) + (k_1718*Ac)/w*(T_L18[n] - T_L17[n])
+ m_dot[n]*Cp*(T_L16[n] - T_L17[n])
+ ka*(Tamb_bot[n] - T_L17[n])*As
)/(M*Cp))
dT_L18_dt = m.Intermediate(((k_1718*Ac)/w*(T_L17[n] - T_L18[n]) + (k_1819*Ac)/w*(T_L19[n] - T_L18[n])
+ m_dot[n]*Cp*(T_L17[n] - T_L18[n])
+ ka*(Tamb_bot[n] - T_L18[n])*As
)/(M*Cp))
dT_L19_dt = m.Intermediate(((k_1819*Ac)/w*(T_L18[n] - T_L19[n]) + (k_1920*Ac)/w*(T_L20[n] - T_L19[n])
+ ka*(Tamb_bot[n] - T_L19[n])*As
)/(M*Cp))
dT_L20_dt = m.Intermediate(((k_1920*Ac)/w*(T_L19[n] - T_L20[n])
+ ka*(Tamb_bot[n] - T_L20[n])*(As+Ac)
)/(M*Cp))
T_L1[n+1] == T_L1[n] + Dt*dT_L1_dt,
T_L2[n+1] == T_L2[n] + Dt*dT_L2_dt,
T_L3[n+1] == T_L3[n] + Dt*dT_L3_dt,
T_L4[n+1] == T_L4[n] + Dt*dT_L4_dt,
T_L5[n+1] == T_L5[n] + Dt*dT_L5_dt,
T_L6[n+1] == T_L6[n] + Dt*dT_L6_dt,
T_L7[n+1] == T_L7[n] + Dt*dT_L7_dt,
T_L8[n+1] == T_L8[n] + Dt*dT_L8_dt,
T_L9[n+1] == T_L9[n] + Dt*dT_L9_dt,
T_L10[n+1] == T_L10[n] + Dt*dT_L10_dt,
T_L11[n+1] == T_L11[n] + Dt*dT_L11_dt,
T_L12[n+1] == T_L12[n] + Dt*dT_L12_dt,
T_L13[n+1] == T_L13[n] + Dt*dT_L13_dt,
T_L14[n+1] == T_L14[n] + Dt*dT_L14_dt,
T_L15[n+1] == T_L15[n] + Dt*dT_L15_dt,
T_L16[n+1] == T_L16[n] + Dt*dT_L16_dt,
T_L17[n+1] == T_L17[n] + Dt*dT_L17_dt,
T_L18[n+1] == T_L18[n] + Dt*dT_L18_dt,
T_L19[n+1] == T_L19[n] + Dt*dT_L19_dt,
T_L20[n+1] == T_L20[n] + Dt*dT_L20_dt
error_top = m.sum((T_top - T_L5[:-1]) ** 2) # Top sensor error
error_mid = m.sum((T_mid - T_L11[:-1]) ** 2) # Mid sensor error
error_bot = m.sum((T_bot - T_L16[:-1]) ** 2) # Bottom sensor error
error_out = m.sum((T_out - T_L18[:-1]) ** 2) # Outlet sensor error
error = error_top + error_mid + error_bot + error_out
m.Minimize(error)
m.solve()
I get the following error:
apm 62.18.232.57_gk_model12 <br><pre> ----------------------------------------------------------------
APMonitor, Version 1.0.1
APMonitor Optimization Suite
--------- APM Model Size ------------
Each time step contains
Objects : 0
Constants : 0
Variables : 0
Intermediates: 0
Connections : 0
Equations : 0
Residuals : 0
@error: Insufficient Model
Model must contain at least 1 variable and 1 equation
STOPPING...
Traceback (most recent call last):
File "xxxx.py", line 214, in <module>
m.solve()
File "C:\Anaconda3\lib\site-packages\gekko\gekko.py", line 2185, in solve
raise Exception(response)
Exception: @error: Insufficient Model
Model must contain at least 1 variable and 1 equation
STOPPING...
I will appreciate if someone could help me here. I have uploaded the dataset file required to run the code at <https://github.com/nagpalh/temporary/blob/main/charge_2.csv>
I checked the .apm file also, it shows the Variables, Intermediates, Connections and Equations noramlly.
</details>
# 答案1
**得分**: 1
使用Gekko的一个很好的应用！尝试使用`IMODE=5`和[有限元上的正交割线][1]，而不是您在积分中使用[Euler方法][2]的`IMODE=3`。目标可能需要修改。4个级别的环境温度是否应与热储罐温度匹配？似乎应该使用4个级别的水的测量温度。4个级别的环境温度应该是模型的输入，用于控制热传递。
[![回归结果][3]][3]
```python
import numpy as np
import pandas as pd
from gekko import GEKKO
import matplotlib.pyplot as plt
data_url = 'https://raw.githubusercontent.com/nagpalh/temporary/main/charge_2.csv'
df = pd.read_csv(data_url)
# 创建GEKKO模型
m = GEKKO()
# 参数
rho = 1000 # 水的密度 [kg/m^3]
r = 0.6 # 储罐的内半径
Ac = np.pi*(r*r) # 储罐的横截面积 [m^2]
Cp = 4186 # 水的比热容 [J/(kg.K)]
H = 1.79 # 储罐的高度 [m]
w = H/20 # 每层的宽度 [m]
M = rho*Ac*w # 每层的质量 [kg]
As = 2*np.pi*r*w # 每层的表面积
Dt = 60 # 采样时间 [秒]
# 输入
m_dot = m.Param() # 质量流量 [kg/s]
T_in = m.Param() # 入口温度 [摄氏度]
Tamb_top = m.Param() # 顶部环境温度 [摄氏度]
Tamb_mid = m.Param() # 中部环境温度 [摄氏度]
Tamb_bot = m.Param() # 底部环境温度 [摄氏度]
m_dot.value = df.mass_flowrate_charge.values
T_in.value = df.Tin_charge.values
Tamb_top.value = df.Tamb_top.values
Tamb_bot.value = df.Tamb_bot.values
Tamb_mid.value = df.Tamb_mid.values
# 观测值
T_out = m.Param() # 出口温度 [摄氏度]
T_top = m.Param() # 顶部传感器温度 [摄氏度]
T_mid = m.Param() # 中部传感器温度 [摄氏度]
T_bot = m.Param() # 底部传感器温度 [摄氏度]
T_out.value = df.Tout_charge.values
T_top.value = df.Ttop_m.values
T_mid.value = df.Tmid_m.values
T_bot.value = df.Tbot_m.values
# 所有层的初始温度
T0 = df.Ttop_m.values
# 初始化数组以存储结果（各层的温度）
T_L1, T_L2, T_L3, T_L4, T_L5, T_L6, T_L7, T_L8, T_L9, T_L10,\
T_L11, T_L12, T_L13, T_L14, T_L15, T_L16, T_L17, T_L18, T_L19, T_L20\
= m.Array(m.Var, 20,value=T0[0])
# 可识别的参数
# k_{i,i+1} = 第i层与第i+1层之间的热导率 [W/(m.K)]. ka - 储罐的热系数 [W/m2]
ka = m.FV(5,lb=0,ub=10000)
p = m.Array(m.FV,19,lb=0.6,ub=10000)
k_12, k_23, k_34, k_45, k_56, k_67, k_78, k_89, k_910, k_1011, k_1112, k_1213,\
k_1314, k_1415, k_1516, k_1617, k_1718, k_1819, k_1920 = p
# 解决方程组
# 每层的能量平衡（层间导热+向环境的热损失+流动层中的流体流动引起的热传递）
m.Equation(M*Cp*T_L1.dt() == (k_12*Ac)/w*(T_L2 - T_L1)
+ ka*(Tamb_top - T_L1)*(As+Ac))
m.Equation(M*Cp*T_L2.dt() == (k_12*Ac)/w*(T_L1 - T_L2) + (k_23*Ac)/w*(T_L3 - T_L2)
+ ka*(Tamb_top - T_L2)*As)
m.Equation(M*Cp*T_L3.dt() == (k_23*Ac)/w*(T_L2 - T_L3) + (k_34*Ac)/w*(T_L4 - T_L3)
+ m_dot*Cp*(T_in - T_L3)
+ ka*(Tamb_top - T_L3)*As)
m.Equation(M*Cp*T_L4.dt() == (k_34*Ac)/w*(T_L3 - T_L4) + (k_45*Ac)/w*(T_L5 - T_L4)
+ m_dot*Cp*(T_L3 - T_L4)
+ ka*(Tamb_top - T_L4)*As)
m.Equation(M*Cp*T_L5.dt() == (k_45*Ac)/w*(T_L4 - T_L5) + (k_56*Ac)/w*(T_L6 - T_L5)
+ m_dot*Cp*(T_L4 - T_L5)
+ ka*(Tamb_top - T_L5))
m.Equation(M*Cp*T_L6.dt() == (k_56*Ac)/w*(T_L5 - T_L6) + (k_67*Ac)/w*(T_L7 - T_L6)
+ m_dot*Cp*(T_L5 - T_L6)
+ ka*(Tamb_top - T_L6)*As)
m.Equation(M*Cp*T_L7.dt() == (k_67*Ac)/w*(T_L6 - T
<details>
<summary>英文:</summary>
Great application of Gekko! Try using `IMODE=5` with [Orthogonal Collocation on Finite Elements][1] instead of `IMODE=3` where you used the [Euler's method][2] for integration. The objective may need to be modified. Should the ambient temperature at the 4 levels match the thermal storage tank temperatures? It seems like it should be using the measured temperatures of the water at the 4 levels. The ambient temperature at the 4 levels should be an input to the model to control the heat transfer.
[![Regression results][3]][3]
```python
import numpy as np
import pandas as pd
from gekko import GEKKO
import matplotlib.pyplot as plt
data_url = 'https://raw.githubusercontent.com/nagpalh/temporary/main/charge_2.csv'
df = pd.read_csv(data_url)
# create GEKKO model
m = GEKKO()
# Parameters
rho = 1000 # density of water [kg/m^3]
r = 0.6 # inner radius of the tank
Ac = np.pi*(r*r) # cross-sectional area of tank [m^2]
Cp = 4186 # specific heat capacity of water [J/(kg.K)]
H = 1.79 # height of the tank [m]
w = H/20 # Width of each layer [m]
M = rho*Ac*w # mass of each layer [kg]
As = 2*np.pi*r*w # surface area of each layer
Dt = 60 # Sample time [sec]
# Inputs
m_dot = m.Param() # mass flow rate [kg/s]
T_in = m.Param() # inlet temperature [deg C]
Tamb_top = m.Param() # ambient temperature at top [deg C]
Tamb_mid = m.Param() # ambient temperature at mid [deg C]
Tamb_bot = m.Param() # ambient temperature at bottom [deg C]
m_dot.value = df.mass_flowrate_charge.values
T_in.value = df.Tin_charge.values
Tamb_top.value = df.Tamb_top.values
Tamb_bot.value = df.Tamb_bot.values
Tamb_mid.value = df.Tamb_mid.values
# Observations
T_out = m.Param() # outlet temperature [deg C]
T_top = m.Param() # top sensor temperature [deg C]
T_mid = m.Param() # mid sensor temperature [deg C]
T_bot = m.Param() # bottom sensor temperature [deg C]
T_out.value = df.Tout_charge.values
T_top.value = df.Ttop_m.values
T_mid.value = df.Tmid_m.values
T_bot.value = df.Tbot_m.values
# Initial temperature for all layers
T0 = df.Ttop_m.values
# Initialize arrays to store results (temperature of layers)
T_L1, T_L2, T_L3, T_L4, T_L5, T_L6, T_L7, T_L8, T_L9, T_L10,\
T_L11, T_L12, T_L13, T_L14, T_L15, T_L16, T_L17, T_L18, T_L19, T_L20\
= m.Array(m.Var, 20,value=T0[0])
# Identifiable parameters
# k_{i,i+1} = thermal conductivity between layer i and i+1 [W/(m.K)]. ka - thermal
# coefficient of tank [W/m2]
ka = m.FV(5,lb=0,ub=10000)
p = m.Array(m.FV,19,lb=0.6,ub=10000)
k_12, k_23, k_34, k_45, k_56, k_67, k_78, k_89, k_910, k_1011, k_1112, k_1213,\
k_1314, k_1415, k_1516, k_1617, k_1718, k_1819, k_1920 = p
# Solve system of equations
# Energy Balance for each layer (conduction between layer + thermal losses
# to environment + heat transfer due to fluid flow in layers)
m.Equation(M*Cp*T_L1.dt() == (k_12*Ac)/w*(T_L2 - T_L1)
+ ka*(Tamb_top - T_L1)*(As+Ac))
m.Equation(M*Cp*T_L2.dt() == (k_12*Ac)/w*(T_L1 - T_L2) + (k_23*Ac)/w*(T_L3 - T_L2)
+ ka*(Tamb_top - T_L2)*As)
m.Equation(M*Cp*T_L3.dt() == (k_23*Ac)/w*(T_L2 - T_L3) + (k_34*Ac)/w*(T_L4 - T_L3)
+ m_dot*Cp*(T_in - T_L3)
+ ka*(Tamb_top - T_L3)*As)
m.Equation(M*Cp*T_L4.dt() == (k_34*Ac)/w*(T_L3 - T_L4) + (k_45*Ac)/w*(T_L5 - T_L4)
+ m_dot*Cp*(T_L3 - T_L4)
+ ka*(Tamb_top - T_L4)*As)
m.Equation(M*Cp*T_L5.dt() == (k_45*Ac)/w*(T_L4 - T_L5) + (k_56*Ac)/w*(T_L6 - T_L5)
+ m_dot*Cp*(T_L4 - T_L5)
+ ka*(Tamb_top - T_L5))
m.Equation(M*Cp*T_L6.dt() == (k_56*Ac)/w*(T_L5 - T_L6) + (k_67*Ac)/w*(T_L7 - T_L6)
+ m_dot*Cp*(T_L5 - T_L6)
+ ka*(Tamb_top - T_L6)*As)
m.Equation(M*Cp*T_L7.dt() == (k_67*Ac)/w*(T_L6 - T_L7) + (k_78*Ac)/w*(T_L8 - T_L7)
+ m_dot*Cp*(T_L6 - T_L7)
+ ka*(Tamb_top - T_L7)*As)
m.Equation(M*Cp*T_L8.dt() == (k_78*Ac)/w*(T_L7 - T_L8) + (k_89*Ac)/w*(T_L9 - T_L8)
+ m_dot*Cp*(T_L7 - T_L8)
+ ka*(Tamb_mid - T_L8)*As)
m.Equation(M*Cp*T_L9.dt() == (k_89*Ac)/w*(T_L8 - T_L9) + (k_910*Ac)/w*(T_L10 - T_L9)
+ m_dot*Cp*(T_L8 - T_L9)
+ ka*(Tamb_mid - T_L9)*As)
m.Equation(M*Cp*T_L10.dt() == (k_910*Ac)/w*(T_L9 - T_L10) + (k_1011*Ac)/w*(T_L11 - T_L10)
+ m_dot*Cp*(T_L9 - T_L10)
+ ka*(Tamb_mid - T_L10)*As)
m.Equation(M*Cp*T_L11.dt() == (k_1011*Ac)/w*(T_L10 - T_L11) + (k_1112*Ac)/w*(T_L12 - T_L11)
+ m_dot*Cp*(T_L10 - T_L11)
+ ka*(Tamb_mid - T_L11)*As)
m.Equation(M*Cp*T_L12.dt() == (k_1112*Ac)/w*(T_L11 - T_L12) + (k_1213*Ac)/w*(T_L13 - T_L12)
+ m_dot*Cp*(T_L11 - T_L12)
+ ka*(Tamb_mid - T_L12)*As)
m.Equation(M*Cp*T_L13.dt() == (k_1213*Ac)/w*(T_L12 - T_L13) + (k_1314*Ac)/w*(T_L14 - T_L13)
+ m_dot*Cp*(T_L12 - T_L13)
+ ka*(Tamb_mid - T_L13)*As)
m.Equation(M*Cp*T_L14.dt() == (k_1314*Ac)/w*(T_L13 - T_L14) + (k_1415*Ac)/w*(T_L15 - T_L14)
+ m_dot*Cp*(T_L13 - T_L14)
+ ka*(Tamb_bot - T_L14)*As)
m.Equation(M*Cp*T_L15.dt() == (k_1415*Ac)/w*(T_L14 - T_L15) + (k_1516*Ac)/w*(T_L16 - T_L15)
+ m_dot*Cp*(T_L14 - T_L15)
+ ka*(Tamb_bot - T_L15)*As)
m.Equation(M*Cp*T_L16.dt() == (k_1516*Ac)/w*(T_L15 - T_L16) + (k_1617*Ac)/w*(T_L17 - T_L16)
+ m_dot*Cp*(T_L15 - T_L16)
+ ka*(Tamb_bot - T_L16)*As)
m.Equation(M*Cp*T_L17.dt() == (k_1617*Ac)/w*(T_L16 - T_L17) + (k_1718*Ac)/w*(T_L18 - T_L17)
+ m_dot*Cp*(T_L16 - T_L17)
+ ka*(Tamb_bot - T_L17)*As)
m.Equation(M*Cp*T_L18.dt() == (k_1718*Ac)/w*(T_L17 - T_L18) + (k_1819*Ac)/w*(T_L19 - T_L18)
+ m_dot*Cp*(T_L17 - T_L18)
+ ka*(Tamb_bot - T_L18)*As)
m.Equation(M*Cp*T_L19.dt() == (k_1819*Ac)/w*(T_L18 - T_L19) + (k_1920*Ac)/w*(T_L20 - T_L19)
+ ka*(Tamb_bot - T_L19)*As)
m.Equation(M*Cp*T_L20.dt() == (k_1920*Ac)/w*(T_L19 - T_L20)
+ ka*(Tamb_bot - T_L20)*(As+Ac))
m.Minimize((T_top - T_L5) ** 2) # Top sensor error
m.Minimize((T_mid - T_L11) ** 2) # Mid sensor error
m.Minimize((T_bot - T_L16) ** 2) # Bottom sensor error
m.Minimize((T_out - T_L18) ** 2) # Outlet sensor error
m.options.IMODE=5
m.time = np.arange(0,len(df))
m.solve()
# open degress of freedom
ka.STATUS = 1
for pi in p:
pi.STATUS = 1
m.options.TIME_SHIFT = 0
m.solve()
# print parameters
print(f'ka: {ka.value[0]}')
for i,pi in enumerate(p):
print(f'p[{i}]: {pi.value[0]}')
# create plots
plt.figure(figsize=(8,5))
plt.subplot(2,2,1)
plt.plot(m.time,T_L1)
plt.plot(m.time,T_L2)
plt.plot(m.time,T_L3)
plt.plot(m.time,T_L4)
plt.plot(m.time,T_L5)
plt.subplot(2,2,2)
plt.plot(m.time,T_L6)
plt.plot(m.time,T_L7)
plt.plot(m.time,T_L8)
plt.plot(m.time,T_L9)
plt.plot(m.time,T_L10)
plt.subplot(2,2,3)
plt.plot(m.time,T_L11)
plt.plot(m.time,T_L12)
plt.plot(m.time,T_L13)
plt.plot(m.time,T_L14)
plt.plot(m.time,T_L15)
plt.subplot(2,2,4)
plt.plot(m.time,T_L16)
plt.plot(m.time,T_L17)
plt.plot(m.time,T_L18)
plt.plot(m.time,T_L19)
plt.plot(m.time,T_L20)
# create plots
plt.figure(figsize=(8,5))
plt.subplot(2,2,1)
plt.plot(m.time,T_L5,'b-',label='T_L5')
plt.plot(m.time,T_top,'rx',label='T_top')
plt.legend()
plt.subplot(2,2,2)
plt.plot(m.time,T_L11,'b-',label='T_L11')
plt.plot(m.time,T_mid,'rx',label='T_mid')
plt.legend()
plt.subplot(2,2,3)
plt.plot(m.time,T_L16,'b-',label='T_L16')
plt.plot(m.time,T_bot,'rx',label='T_bot')
plt.legend()
plt.subplot(2,2,4)
plt.plot(m.time,T_L18,'b-',label='T_L18')
plt.plot(m.time,T_out,'rx',label='T_out')
plt.legend()
plt.show()
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库，让每个人都能够通过互相帮助和分享经验来进步。
评论