这几天在学习WOFOST模型的过程中遇到了一个问题:对于根据时期变化的变量该怎么设置才可以成功进行 SOBOL 敏感性分析和 SUBPLEX 参数优化。苦于国内关于 PCSE 的资料甚少,大部分都是直接照搬官方文档,而官方文档并没有这方面的示例,困扰了挺久后问了一位大佬让我有了新思路,多亏了他才有我这篇文章。
由于本人资质尚浅,某些想法和文章中的话语难免有些幼稚,这个望各位不啬指正。
先上大佬的文章,他是利用simlab软件进行敏感性分析,然后用图片中的方法进行模型参数替换。
但是我想用代码的方式实现 SOBOL 方法进行敏感性分析。还有就是本人想用面向对象的写法实现模型,尽量的精简代码,避免遇到不同的作物、不同的情况又得频繁的修改代码。所以在这里我详细地为大家介绍我的方法:
首先对于根据时期变化的变量我们可以分成多个变量
以 AMAXTB 为例,下图为官方示例的小麦作物参数文件
可以看见这有四个数值,分别对应发育期为0.00、1.00、1.30、2.00的 AMAXTB 值。那我们在设置参数范围时,可以这样:
分成四个变量,后面的数字代表发育阶段。每个变量对应四个范围(在这里我是四个[20,70])。
注意这个变量后的数字要和作物参数文件上的发育阶段表示方法一致,不然我后面修改的代码会无法成功匹配和替换。
设置完成后读取这个参数范围yaml文件,在敏感性分析代码(我将敏感性分析写成了一个单独的模块,各位根据自己的情况修改)中添加一个判断:
def Sensitivity(self,sampledata: str):
# 使用 Saltelli 采样器采集生成参数集
calc_second_order = True
nsamples = ((65*(self.problem['num_vars'])//500)+1)*500
print("样本数为 {}".format(nsamples))
paramsets = saltelli.sample(self.problem, nsamples, calc_second_order=calc_second_order)
print("We are going to do %s simulations" % len(paramsets))
nsets = len(paramsets)
printProgressBar(0, nsets, prefix = 'Progress:', suffix = 'Complete', length = 50)
# 循环遍历参数值的样本,并为每个样本运行 WOFOST
target_results = []
for i, paramset in enumerate(paramsets):
self.params.clear_override() # 清除任何现有的覆盖
for name, value in zip(self.problem["names"], paramset):
pername = re.findall(r"\d+\.\d+|[A-Z]+",name)
if pername[0] in ["SLATB","KDIFTB","EFFTB","AMAXTB","TMPFTB"]:
index = self.params[pername[0]].index(float(pername[1]))+1
self.params[pername[0]][index] = value
else:
self.params.set_override(name, value) # 更改每个 WOFOST 模拟的模型参数
wofost = Wofost72_PP(self.params, self.wdp, self.agro) # 导入模型
wofost.run_till_terminate() # 运行模型
r = wofost.get_summary_output() # 返回已在模拟期间存储的汇总变量。
target_result = r[0][self.target_variable]
if target_result is None:
print("Target variable is not available in summary output!")
break
target_results.append(target_result)
printProgressBar(i+1, nsets, prefix = 'Progress:', suffix = 'Complete', length = 50)
target_results = np.array(target_results)
# 计算灵敏度指数,使用 Sobol 方法来评估模型对参数变化的响应。
Si = sobol.analyze(self.problem, target_results, calc_second_order=calc_second_order)
这句代码是利用正则(re库)对变量进行分割变成字符串和数字(包括小数)。这样做的目的是分割出来的字符串模型是可以正常识别的,而数字对应变量的发育阶段,之后可以依据这个数字进行值替换。
pername = re.findall(r"\d+\.\d+|[A-Z]+",name)
这句代码是进行判断,是否是根据时期变化的变量。列表中的变量是一些根据时期变化的变量,如果有其他可以自己增加。
if pername[0] in ["SLATB","KDIFTB","EFFTB","AMAXTB","TMPFTB"]:
这句代码是根据发育阶段得到该进行值替换的索引
index = self.params[pername[0]].index(float(pername[1]))+1
最后就可以成功分情况进行值替换啦。
self.params[pername[0]][index] = value
最后再次感谢大佬的启发!希望其他在学习WOFOST的小伙伴看到我这篇文章可以有所收获!