import numpy as np
import matplotlib.pyplot as plt
from sklearn.metrics import mean_squared_error

基本参数

initial_infected_leaves = 1 # 初始感染叶片数
initial_susceptible_leaves = 1000 # 初始健康叶片数
days = 4 # 模拟时间段（天数）

计算每日感染率的函数

def calculate_daily_infection_rate(DPi, RAi, WIi, DPi_sqrt, DTi_prime):
ln_R = (-0.07988 * DPi) + (0.09983 * RAi) + (0.6276 * np.log(WIi)) + (0.7448 * DPi_sqrt) + (0.06967 * DTi_prime * DPi) + 0.8616
return np.exp(ln_R)

显症率函数

def calculate_symptom_rate(latent_leaves, symptom_rate_param):
symptom_leaves = latent_leaves * symptom_rate_param
return symptom_leaves

病斑扩展率函数

def calculate_disease_expansion_rate(symptom_leaves, expansion_rate_param):
expanded_leaves = symptom_leaves * expansion_rate_param
return expanded_leaves

报废率函数

def calculate_discard_rate(disease_leaves, discard_rate_param):
discarded_leaves = disease_leaves * discard_rate_param
return discarded_leaves

示例输入（应使用真实数据）

DPi = 3 # 露时（小时）
RAi = 0.1 # 降雨量（毫米）
WIi = 1 # 风速（米/秒）
DPi_sqrt = np.sqrt(DPi)
DTi_prime = 1 # 温度生长当量

计算每日感染率

daily_infection_rate = calculate_daily_infection_rate(DPi, RAi, WIi, DPi_sqrt, DTi_prime)
print(f"每日感染率: {daily_infection_rate}")

初始化数组以存储每天的值

susceptible_leaves = np.zeros(days)
infected_leaves = np.zeros(days)
new_infections = np.zeros(days)
symptom_leaves = np.zeros(days)
expanded_leaves = np.zeros(days)
discarded_leaves = np.zeros(days)

设置初始值

susceptible_leaves[0] = initial_susceptible_leaves
infected_leaves[0] = initial_infected_leaves

显症率、扩展率和报废率的假设参数

symptom_rate_param = 0.3
expansion_rate_param = 0.4
discard_rate_param = 0.2

模拟循环

for day in range(1, days):
# 计算新感染叶片数
new_infections[day] = daily_infection_rate * infected_leaves[day - 1] * (susceptible_leaves[day - 1] / (susceptible_leaves[day - 1] + infected_leaves[day - 1]))

# 计算显症叶片数
symptom_leaves[day] = calculate_symptom_rate(new_infections[day], symptom_rate_param)

# 计算病斑扩展叶片数
expanded_leaves[day] = calculate_disease_expansion_rate(symptom_leaves[day], expansion_rate_param)

# 计算报废叶片数
discarded_leaves[day] = calculate_discard_rate(expanded_leaves[day], discard_rate_param)

# 更新健康叶片数和感染叶片数
susceptible_leaves[day] = susceptible_leaves[day - 1] - new_infections[day]
infected_leaves[day] = infected_leaves[day - 1] + new_infections[day]

绘制结果

plt.figure(figsize=(10, 6))
plt.plot(range(days), susceptible_leaves, label="susceptible_leaves")
plt.plot(range(days), infected_leaves, label="infected_leaves")
plt.plot(range(days), symptom_leaves, label="symptom_leaves")
plt.plot(range(days), expanded_leaves, label="expanded_leaves")
plt.plot(range(days), discarded_leaves, label="discarded_leaves")
plt.xlabel("Days")
plt.ylabel("Leaf number")
plt.legend()
plt.title("Wheat stripe rust - Comprehensive Model")
plt.show()

假设我们有一些观测数据用于验证

observed_infected_leaves = np.array([1, 20, 45, 400])

绘制预测值与观测值的对比

plt.figure(figsize=(10, 6))
plt.plot(range(days), infected_leaves, label="Predicted")
plt.plot(range(len(observed_infected_leaves)), observed_infected_leaves, 'o', label="Observed")
plt.xlabel("Days")
plt.ylabel("Infected leaf number")
plt.legend()
plt.title("Model Test")
plt.show()

计算误差指标，如均方误差

mse = mean_squared_error(observed_infected_leaves, infected_leaves[:len(observed_infected_leaves)])
print(f"MSE: {mse}")