综合能源系统优化建模：P2G+CCS+CHP耦合的Matlab实现

最近在做综合能源系统方向的项目核心内容就是把碳捕集系统CCS、电转气P2G和热电联产CHP这三块耦合到同一个优化模型里用Matlab做建模、求解和仿真分析。这个方向这两年热度确实高因为它同时踩中了两个痛点一是高比例新能源并网后的弃风弃光问题二是双碳目标下的碳排放约束问题。这篇博文我打算把一个完整的P2GCCSCHP综合能源系统建模优化流程拆开讲清楚——从物理过程怎么翻译成数学公式到优化模型怎么搭再到Matlab代码怎么组织、求解器怎么选最后把我实际调参仿真时踩过的坑也一并列出来。适合正在做综合能源系统优化、或者刚接触碳捕集/电转气建模的研究生和工程师参考目标是你看完能直接照着搭出一个可运行的框架。1. 为什么要把P2G、碳捕集和热电联产放在一个模型里1.1 单看CHP是稳定器放进高比例新能源系统就成了麻烦源传统热电联产机组最大的特点是以热定电冬季供热期为了满足热负荷机组必须维持较高的电出力而且供热出力越大电出力下限越高。这个物理特性在纯火电系统里问题不大但一旦系统里接入大规模风电、光伏麻烦就来了——热负荷高峰往往和风电大发时段重叠CHP被迫满发电网消纳空间被挤占风电只能白白弃掉。我在实际仿真里见过供热季弃风率超过20%的算例而且越是供暖季弃风越严重这就是典型的热电耦合困局。解决思路无非两条一是给CHP减负让它在风电大发时少发电二是给系统加荷把多余的电力变成其他形式的能量存起来或利用掉。前者对应的是改造热电联产运行方式比如加装储热罐、电锅炉后者对应的就是电转气技术。1.2 P2G和CCS各自解决了什么耦合起来又解决了什么电转气Power-to-Gas简单说就是把电能转化为燃气通常分两步第一步电解水制氢第二步氢气与CO2甲烷化生成天然气。注意第二步需要消耗CO2而这个CO2从哪来来如果单独买成本很高如果和碳捕集系统配合CCS正好从烟气里捕集CO2直接作为P2G原料——这就形成了一个很有意思的闭环CHP燃烧天然气发电产热产生含CO2烟气CCS从烟气中捕集CO2一部分外送封存或利用一部分供给P2GP2G利用风电低谷期的多余电能电解水制氢再和CO2甲烷化产出天然气回补气网天然气又可以供给CHP使用。这个循环的好处不止一个。从消纳端看P2G相当于一个可调节的大功率负荷风电大发时多开风电不足时少开直接抬高了系统的灵活调节能力。从碳端看CCS把本来要排到大气里的CO2固定下来P2G又把它变成了燃料碳排放强度显著下降。从热端看CCS的溶剂再生过程需要消耗大量蒸汽而这部分蒸汽通常从CHP的抽汽中获取相当于给CHP的对外供热能力加了约束反而缓解了以热定电——机组不能把全部抽汽都拿去供热了电出力下限跟着松动给风电腾出了更多空间。1.3 三设备耦合建模的本质所以把P2G、CCS和CHP放在同一个模型里本质上是把三条能量流电能流、热能流、天然气流和一条碳流CO2捕集、利用、排放统一建模。电力平衡约束里P2G是新增电负荷热力平衡约束里CCS是新增热负荷储热罐可以削峰填谷燃气平衡里P2G产气是气源碳约束里CCS捕集量和P2G消耗量是衔接变量。四个平衡关系通过耦合变量互相牵制这才是这个优化问题的核心难点也是Matlab建模中最需要仔细处理的地方。2. 系统架构与设备物理模型能量流怎么走约束从哪来2.1 系统拓扑和能量流关系我在模型里设置的综合能源系统结构大概是这样的风电接入电网CHP同时接入电网和热网燃气锅炉作为热网备用热源电锅炉可以直接用电供热P2G从电网取电产气并接入天然气网CCS从CHP烟气中捕集CO2部分供给P2G储热罐挂在热网上用于时段转移。电、热、气三种能源通过各自的母线平衡耦合点就在CHP机组、P2G和CCS这三个设备上。做模型第一步就是把所有设备的输入输出关系列成一张表CHP输入天然气输出电能和热能P2G输入电能和CO2输出天然气CCS输入热能和烟气输出CO2流储热罐输入/输出热能电锅炉输入电能输出热能。这张表就是后面写约束方程的基础建议新手第一步先画这张表而不是急着写代码。我见过太多人上来就写Yalmip代码结果能量流方向搞反了后面调试极其痛苦。2.2 CHP机组建模可行域而不是单点热电联产机组的关键在于电出力和热出力不是独立的而是有一个耦合可行域。对于最常见的抽凝式机组可行域可以用一组线性不等式来逼近描述。我用的简化模型是这样的% CHP电出力与热出力耦合约束抽凝式 % P_e_min c_v*P_h P_e P_e_max - k_v*P_h % 其中 c_v 是低压供热工况下的电出力下限斜率k_v 是电出力上限斜率 P_e_min 20; % MW 纯凝工况最低电出力 P_e_max 200; % MW 纯凝工况最大电出力 c_v 0.15; % 供热对电出力下限的影响系数 k_v 0.25; % 供热对电出力上限的影响系数 P_h_max 120; % MW 最大热出力这里每个参数都有明确的物理意义当热出力增加时电出力上限会下降因为抽汽多了做功的蒸汽少了电出力下限会上升为了维持供热机组不能降得太低。这两条线的斜率k_v和c_v是CHP机型的关键参数来自机组的实际运行数据仿真时不能随意拍脑袋最好参考同类型机组的铭牌参数或文献数据。CHP的燃料消耗量我用线性化模型即燃料成本是电出力和热出力的线性组合% F_chp(t) b0 b1*P_e(t) b2*P_h(t) b0 40; % 空载燃料成本 (万元/h) b1 0.25; % 电出力燃料成本系数 (万元/MWh) b2 0.12; % 热出力燃料成本系数 (万元/MWh)注意这个线性化只在机组运行区间内有效不能外推到可行域之外。如果你的项目需要更高精度可以把二次成本函数做分段线性化后文第6章会专门讲这个问题。2.3 P2G模型电解槽效率与甲烷化的CO2需求P2G我建模为两段电解水制氢电→氢和甲烷化氢CO2→甲烷水。整体电转气效率指的是电→甲烷的效率即G_p2g(t) eta_p2g * P_p2g(t) * 3.6 / LHV_gas; % eta_p2g 综合效率约 0.45~0.60电转甲烷 % LHV_gas 为天然气低位热值约为 36 MJ/m3 % P_p2g 单位 MWG_p2g 单位 km3/h千立方米/小时这个式子为什么是3.6因为1 MWh电能量等于3.6 GJ需要除以天然气的低位热值才能换算成体积。很多新手在这里忘记单位换算算出来的产气量差了好几个数量级后面所有平衡约束全乱套。甲烷化过程的CO2需求量按化学计量比来算每生产1立方米甲烷大约需要1.8公斤CO2。写成约束就是% 每产 1 km3 甲烷需要的CO2质量吨 CO2_demand(t) 1.8 * G_p2g(t); % 单位吨/h % 注意1 km3 1000 m31.8 kg/m3 1.8 t/km3这个计算看起来简单但它在模型里起着关键的耦合作用CCS捕集量不仅要覆盖系统自身的碳排放控制目标还要保证P2G的CO2原料供应充足。如果CCS容量太小P2G就没有原料只能额外购买CO2经济性直线下降。2.4 CCS模型捕集率、能耗和热负荷耦合碳捕集系统我按燃烧后胺法捕集来建模核心参数是捕集效率和单位能耗。捕集到的CO2量和烟气中CO2总量之间的关系是P_capture(t) alpha_ccs * (E_chp_emission(t) - E_vent(t)); % alpha_ccs 捕集率一般取 0.85~0.95 % E_chp_emission 是CHP单位时段的原始CO2排放量吨/h % P_capture 是实际捕集到的CO2量吨/hCCS系统最大的特点是再生塔需要大量蒸汽来加热富液解析CO2这部分蒸汽通常从CHP的中低压抽汽中获取。我用的模型是H_ccs(t) theta_ccs * P_capture(t); % theta_ccs 为单位捕集的热耗系数典型值 0.2~0.4 MWh/tCO2这里最关键的地方在于CCS消耗的热能会从CHP对外供热中扣除。我在热平衡约束里写的是P_h(t) - H_ccs(t)才是实际对外供热量这意味着CCS捕集量越大机组能供给热网的热量越少间接拉低了CHP必须维持在较高出力的压力。这个机理在实际论文里经常被描述为碳捕集削弱热电耦合是碳捕集系统对新能源消纳起到积极作用的核心所在。2.5 储能系统模型储热罐的模型相对简单就是能量守恒加上容量约束S_op(t) S_op(t-1) (H_charge(t) - H_discharge(t)) * dt; 0 S_op(t) S_max; H_charge(t) H_charge_max; H_discharge(t) H_discharge_max;我要特别提醒一个细节储热罐的初末状态约束。实际工程中调度是周期性的通常要求一个调度周期结束时储热量回到初始值否则模型会把储热罐当成免费的无限能量库导致结果失真。我在代码里加了S_op(1) S_op(T)作为周期约束这个约束必须在建模阶段就加上。3. 优化问题构建目标函数、等式不等式约束与求解器取舍3.1 目标函数为什么我选择运行成本最小化综合能源系统优化的目标函数一般有几个选择最小化总运行成本、最小化碳排放、最小化弃风量或者做多目标加权。我在这个项目里选的是运行成本最小化同时把碳排放通过碳交易机制折算进成本原因有三第一成本最小化是最接近实际运营者的决策逻辑工程落地性最强。第二碳排放通过碳价进入目标函数后不需要额外处理多目标权重模型还是单目标优化求解简单且结果好解释。第三仿真想回答的核心问题是CCSP2G在经济上是否划算只有把成本作为结果指标才能对比不同技术组合的可行性。目标函数的具体构成% 总成本 购气成本 购电成本 碳交易成本 设备运维成本 Objective sum(c_gas * F_chp c_gas * F_gb) ... % 燃气成本 sum(c_grid * F_net) ... % 购电成本 sum(c_carb * (E_total - E_allowance)) ... % 碳交易成本 sum(c_om_chp * P_e c_om_p2g * P_p2g c_om_ccs * P_capture); % 运维成本注意碳交易成本一项里E_total是系统总排放E_allowance是免费配额两者之差乘以碳价就是实际碳成本。这个模型里如果CCS捕集量足够大E_total - E_allowance可以变成负数也就是系统通过卖碳配额赚钱这在实际运行中是真实存在的激励机制。3.2 约束条件清单我把模型里的约束分成五类写代码时这样分类有助于排查问题首先是平衡约束。电功率平衡、热功率平衡、天然气平衡、CO2平衡四条属于等式约束直接反映能量守恒。电功率平衡的典型写法是P_e(t) P_wind(t) F_net(t) E_load(t) P_eb(t) P_p2g(t) P_hes(t); % P_hes 是电储能充电功率如果模型里有电储能的话热功率平衡则是P_h(t) H_gb(t) H_discharge(t) H_load(t) H_charge(t) H_ccs(t);第二是设备运行约束包括CHP可行域、锅炉出力范围、P2G和CCS的容量限制。第三是爬坡约束CHP电出力相邻时段变化率一般限制在每小时5%~10%额定容量以内。第四是储能动态约束。第五是碳排放约束包括总碳排放不能超过上限、CCS捕集量受限于烟气中CO2总量等。3.3 求解器取舍MILP是默认选择确定目标函数和约束之后接下来要决定模型用什么样的求解器。我的建议是除非有强烈的非线性环节必须保留否则一律把模型线性化用混合整数线性规划MILP求解。原因很实在综合能源系统调度本质上是多时段耦合决策问题典型周期24小时或168小时变量规模几千到几万个非线性求解器如fmincon、ipopt很容易陷入局部最优而且求解时间不可控。而线性/MILP模型配合Gurobi或CPLEX在求解质量和速度上都远优于通用非线性求解器。我在Matlab里的标准配置是Yalmip做建模层Gurobi做底层求解器ICPlex也可以效果类似。% Yalmip设置求解器 ops sdpsettings(solver, gurobi, verbose, 2); solution optimize(Constraints, Objective, ops);如果你的电脑没有装Gurobi退而求其次可以用MATLAB自带的intlinprog但求解大规模问题时会明显慢一截。我建议课题组有条件的话用学术版Gurobi免费而且性能极佳安装也就是注册一个academic license的事。4. Matlab代码实现从参数表到可运行的优化主程序4.1 代码结构怎么组织我做这类项目的Matlab代码通常分四个文件避免全部逻辑堆在一个脚本里init_params.m定义所有设备参数、预测数据风电出力曲线、负荷曲线、价格参数输出一个结构体数组pbuild_model.m接收参数p声明Yalmip变量、写约束和目标函数返回Constraints和Objectivesolve_case.m调用optimize求解做后处理把结果整理成表格和曲线plot_results.m画图包括电平衡曲线、热平衡曲线、弃风情况、碳排放分布。我目前用的是Matlab 2024bYalmip需要安装最新版推荐从GitHub拉R2024版本。需要提醒一下Matlab版本和Yalmip的兼容性有时候会出问题特别是新版Matlab里Yalmip的路径配置如果没弄好会报一些奇怪的错误。还有就是环境变量的配置问题比如Win10下装了多个Matlab版本启动时容易遇到mathworks licensing error 8之类的许可证报错这种基本都是许可证文件路径或者MLM_LICENSE_FILE环境变量没配好导致的和代码本身没关系遇到先检查环境别急着怀疑模型。4.2 主程序搭建变量声明和核心约束下面这段代码是build_model.m的核心部分我把变量声明和关键约束写出来方便你直接对照着搭建自己的模型function [Constraints, Objective, variables] build_model(p) T p.T; % 调度时段数例如24 %% 决策变量 % CHP P_e sdpvar(1, T); % 电出力 MW P_h sdpvar(1, T); % 热出力 MW F_chp sdpvar(1, T); % 燃料量折合天然气 k m3/h % P2G P_p2g sdpvar(1, T); % 电耗 MW G_p2g sdpvar(1, T); % 产天然气 k m3/h % CCS P_cap sdpvar(1, T); % 实际捕集CO2量 t/h H_ccs sdpvar(1, T); % 热耗 MW % 电网 F_net sdpvar(1, T); % 网购电 MW % 锅炉、储能 H_gb sdpvar(1, T); % 燃气锅炉热出力 MW H_ch sdpvar(1, T); % 储热充电功率 MW H_dis sdpvar(1, T); % 储热放电功率 MW S_op sdpvar(1, T); % 储热罐蓄电量 MWh %% 约束 Constraints []; % 1) CHP可行域 Constraints [Constraints, p.P_e_min p.c_v * P_h P_e p.P_e_max - p.k_v * P_h]; Constraints [Constraints, 0 P_h p.P_h_max]; % CHP燃料线性化成本 Constraints [Constraints, F_chp p.b0 p.b1 * P_e p.b2 * P_h]; % 2) P2G Constraints [Constraints, G_p2g p.eta_p2g * P_p2g * 3.6 / p.LHV_gas]; Constraints [Constraints, 0 P_p2g p.P_p2g_max]; % 3) CCS Constraints [Constraints, P_cap p.alpha_ccs * p.CO2_coef * F_chp]; % 捕集量不能超过排放量这里 CO2_coef 为单位燃气燃烧排放系数 Constraints [Constraints, P_cap p.CO2_coef * F_chp]; % CCS热耗 Constraints [Constraints, H_ccs p.theta_ccs * P_cap]; % 4) 电功率平衡 Constraints [Constraints, P_e p.P_wind F_net p.E_load P_p2g]; % 5) 热功率平衡 Constraints [Constraints, P_h H_gb H_dis p.H_load H_ch H_ccs]; % 6) 储热 Constraints [Constraints, S_op(1) p.S_init]; for t 2:T Constraints [Constraints, S_op(t) S_op(t-1) (H_ch(t) - H_dis(t)) * p.dt]; end Constraints [Constraints, S_op(T) p.S_init]; % 周期性 Constraints [Constraints, 0 S_op p.S_max]; Constraints [Constraints, 0 H_ch p.H_ch_max, 0 H_dis p.H_dis_max]; % 7) 锅炉 Constraints [Constraints, 0 H_gb p.H_gb_max]; %% 目标函数 Objective sum(p.c_gas * F_chp p.c_gas * F_gb) ... sum(p.c_grid * F_net) ... sum(p.c_carb * (p.CO2_coef * F_chp - p.alpha_ccs * P_cap - p.E_allowance)) ... sum(p.c_om_chp * P_e p.c_om_p2g * P_p2g); %% 返回变量用于后处理 variables struct(P_e, P_e, P_h, P_h, P_p2g, P_p2g, ... P_cap, P_cap, F_net, F_net, S_op, S_op); end写这个模型有个小经验Yalmip里sdpvar声明向量变量时下标从1到T循环写约束时注意Matlab的索引习惯用optimize求解完以后直接用value(P_e)取值别再用value(P_e)这种字符串方式新版Yalmip对字符串变量名的支持已经移除了旧代码经常在这里报错。4.3 求解和后处理求解部分很简单就是一个optimize调用。但后处理我认为比求解本身更花心思。我通常把结果整理成一张表每个时段的风电消纳量、CHP电出力、P2G电耗、CCS捕集量、储热罐SOC、购电量、碳排放。然后画6张图电平衡堆叠面积图、热平衡堆叠面积图、风电消纳曲线、碳排放柱状图、储热SOC曲线和成本构成饼图。这些图是论文里最常用的结果图也是你判断模型是否合理的直接依据——比如某时段储热SOC曲线震荡剧烈的多半是充放策略出了问题要回头查约束。5. 典型日仿真三种运行模式下的成本、碳排放与弃风率对比5.1 场景设置为了回答CCS和P2G到底带来了什么我给模型设计了三个运行模式做对照实验。模式A是基准系统没有CCS也没有P2G只有CHP、锅炉、储热和风电模式B在A基础上加了CCS模式C是完整系统CCS和P2G都投入运行。风电出力曲线和电、热负荷曲线都采用同一个典型冬季日数据风电装机150MW电负荷峰值180MW热负荷峰值120MWCHP容量200MW。这里有个参数设置上的小细节风电出力和热负荷的匹配关系在冬季很微妙。白天热负荷高、风电一般夜间热负荷依然高因为是供热季风电却常常达到峰值这就是弃风发生的典型窗口。如果你的仿真结果里完全没有弃风时段说明风电装机容量太小系统不具备分析价值建议先调大风电渗透率。5.2 仿真结果与对比我在一组典型参数下跑出来的结果大致如下具体数值随参数略有浮动这里给出的是趋势性结果指标模式A无CCS无P2G模式B含CCS模式C含CCSP2G总运行成本万元/日82.596.391.8碳排放量吨/日486289214弃风率%18.612.25.4风电上网电量MWh/日158017101840P2G耗电量MWh/日00310CCS捕集量吨/日0176198从结果能看出几个有意思的规律。第一模式B相比模式A碳排放下降明显从486降到289吨但总成本反而上升了13.8万。这是因为CCS运行需要消耗大量蒸汽机组燃料成本增加而且碳交易节省的钱在低碳价下还不足以覆盖CCS的额外能耗成本。这说明加装CCS是有成本的不能只想着环保不看钱。第二模式C在模式B的基础上增加了P2G碳排放又降了75吨弃风率从12.2%降到5.4%成本反而比模式B降了4.5万。原因在于P2G把原本要弃掉的风电转化为天然气这部分燃气是有价值的抵消了P2G设备本身的运维成本。也就是说弃风严重时P2G不仅环保而且经济上也能回血。第三仔细观察模式C的调度策略会发现P2G基本都在凌晨风电高峰时段满功率运行白天风电不足时功率降低甚至停机CCS则在夜间CHP高负荷时段捕集量最大白天捕集量减少。这种时段差异是优化模型自己学出来的体现了多能互补的价值。5.3 碳价敏感性分析我顺手做了个碳价敏感性测试把碳价从每吨30元调到300元模式C的优势越来越明显。碳价低的时候模式B和C总成本高于模式A但当碳价超过每吨150元时模式C的总成本开始低于模式A也就是低碳技术在经济上self-validating。这个结论在项目汇报里很有用可以论证碳交易机制是P2GCCS组合可行性的关键支撑。6. 调试阶段最常犯的五个错误与排查经验6.1 单位换算错误最隐蔽也最致命这个方向建模最坑的就是单位。MW、MWh、GJ、km3、吨、公斤稍微一个换算错结果就离谱。我自己最惨的一次是把LHV_gas的单位搞错导致P2G产气量算大了1000倍天然气平衡完全失衡模型直接不可行。排查了整整一天才意识到是低位热值36MJ/m3和36GJ/km3之间的换算出了问题。我现在的做法是所有单位换算在参数初始化脚本里集中处理并用注释标清楚每一步的单位。比如% 1 km3 天然气 1000 m3 % LHV_gas 36 MJ/m3 36e3 MJ/km3 36 GJ/km3 % 1 MWh 3.6 GJ然后专门写一个unit_check函数测试几个已知数值比如1MWh电通过效率0.5的P2G应该产出多少天然气先验证再跑优化绝不跳过。6.2 初末状态约束忘了加储能变成能量印钞机如果不加S_op(1) S_op(T)这个周期约束优化结果里储热罐会在最后一个时段把全部热量放光因为储热罐剩余热量不计入目标函数白白放弃不利用就是浪费。加了周期约束后储热罐的初末蓄热量保持一致结果才符合实际日调度逻辑。这个问题几乎每个新手都会踩我在第2.5节已经强调过这里再次列出是因为太常发生了。6.3 线性化分段的精度问题CHP的燃料成本如果直接用二次函数F a*P^2 b*P c优化问题就变成二次规划QP虽然Gurobi也能解二次问题但连续二次模型对拟合数据的依赖很强而且二次函数在整个区间单调性可能与实际不符。我的方案是分段线性化把电出力区间切成4~6段每段用不同的线性斜率和截距表示燃料消耗再引入二进制变量选择区间。这样模型是MILP求解稳定精度也够用。分段线性化需要引入0-1整数变量会增加求解时间但24小时场景下影响不大。如果你跑168小时一周的模型觉得慢优先检查是不是整数变量太多而不是急着换求解器。6.4 热平衡方程里漏掉CCS热耗模式B刚搭好时我一度发现结果很奇怪加了CCS之后CHP热出力并没有明显增加碳排放也没有预期的降幅。后来才发现热平衡方程里只写了P_h H_load H_ch漏掉了 H_ccs这一项。CCS消耗的热能没有进热平衡约束等于模型认为CCS不耗能捕集CO2是免费的——这当然是错的。这个错很好查只需要对比加装CCS前后CHP总热出力的变化如果几乎不变那一定有约束漏写了。6.5 风水曲线和负荷曲线时间对齐问题调数据时还有个容易忽略的坑风电出力曲线通常是15分钟或1小时粒度热负荷曲线可能是15分钟粒度如果直接把两个不同采样率的数据塞进同一个24×1的向量里维度对不上Yalmip会报维度不一致错误。遇到这种问题先统一时间粒度我通常全部重采样到1小时再用线性插值补齐。还要注意时区问题特别是跨夏令时的预测数据虽然国内没有夏令时但如果用的是国外公开数据集就要仔细核对时间轴。7. 几个我认为值得多花时间的扩展方向模型做到能跑通、出结果只是第一步。我自己的项目里有几个方向明显值得继续投入精力这里也一并分享给正在规划后续工作的同行一点参考。第一个是考虑更精细的P2G模型。目前我把电解槽和甲烷化看成一个整体效率恒定实际上电解槽在部分负荷下的效率会下降启动和停机也有动态响应约束。如果要研究P2G参与调频或备用市场的价值恒效率模型会给出偏乐观的结果需要改成变效率曲线或分段线性模型。第二个是引入需求响应。电、热负荷可以设置一部分可中断/可转移负荷把它们作为灵活性资源参与调度。这样系统的弹性会显著提高P2G和储热罐的容量需求可能下降经济性结论会有变化。第三个是考虑碳排放的隐含成本。我这里用的是碳交易价格属于市场机制。如果项目面向的是碳约束政策场景可以把碳排放上限作为硬约束然后对比不同排放目标下的系统成本和弃风率输出一条碳排放-成本帕累托曲线这个图在论文和项目报告中都非常好用。我在实际操作中最深的体会是这套综合能源系统建模的核心价值不是解出一个最优解而是让你真正理解电、热、气、碳四条能量流之间的耦合逻辑。第一次跑通模式C的时候看着储热罐在夜间蓄热、白天放热P2G在风电高峰启动CCS捕集量跟着CHP出力走你会真切感受到系统优化这四个字的含义。MatlabYalmipGurobi这套技术栈足够支撑你在这个方向上做很深的探索关键是每一步都搞清楚物理过程再翻译成数学语言。做仿真研究模型对了结论才有意义。