import pandas as pd
import numpy as np
from openpyxl import load_workbook
# 参数设置
input_excel_file = 'output_data.xlsx'
output_excel_file = '附件6-问题二答案表.xlsx'
rho = 1.225 # 空气密度 (kg/m^3)
A = 5026 # 截面积 (m^2)
num_turbines_per_sheet = 100 # 每个工作表风机数量
num_time_steps = 2000
Cp=0.593
# 工作表名称与对应的数据源
sheet_config = {
'主轴扭矩': 'WF_1',
'塔架推力': 'WF_2'
}
# 构建时间列
time_column = np.arange(1, num_time_steps + 1).reshape(-1, 1) # (2000, 1)
# 用于保存结果
results = {
'主轴扭矩': [],
'塔架推力': []
}
# 遍历每个工作表及其对应的数据源
for sheet_name, source_sheet in sheet_config.items():
print(f"正在处理工作表: {sheet_name}")
# 一次性读取整个数据源工作表
df = pd.read_excel(input_excel_file, sheet_name=source_sheet, header=None)
data = df.values # 转为 numpy 数组
# 处理每台风机
for turbine_idx in range(num_turbines_per_sheet):
start_row = turbine_idx * num_time_steps
end_row = start_row + num_time_steps
# 提取当前风机的数据块
turbine_data = data[start_row:end_row, :]
# 提取第二列为有功功率 P,第三列为等效风速 v
p = turbine_data[:, 1]
v = turbine_data[:, 2]
# 计算应力或扭矩
valid = v != 0
stress_or_torque = np.zeros_like(v, dtype=float)
stress_or_torque[valid] = (0.5 * rho * v[valid]**3 * A*Cp - p[valid]) / v[valid]
# 存入结果列表
results[sheet_name].append(stress_or_torque)
# 写入 Excel 文件
with pd.ExcelWriter(output_excel_file, engine='openpyxl', mode='a', if_sheet_exists='replace') as writer:
for sheet_name in sheet_config:
# 将所有风机数据拼接成 (2000, 100)
turbines_data = np.column_stack(results[sheet_name]) # 合并为 100 列
# 添加时间列
final_df = pd.DataFrame(np.hstack([time_column, turbines_data]),
columns=['时间'] + [f'风机{i+1}' for i in range(num_turbines_per_sheet)])
# 写入对应的工作表
final_df.to_excel(writer, sheet_name=sheet_name, index=False)
print(f"应力/扭矩数据已成功写入文件:{output_excel_file}")
根据此代码中计算应力/扭矩的模型,对以下代码进行修改
import numpy as np
import pandas as pd
from pykalman import KalmanFilter
from scipy.interpolate import interp1d
from deap import base, creator, tools, algorithms
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
plt.rcParams['font.sans-serif'] = ['SimHei'] # 中文显示
plt.rcParams['axes.unicode_minus'] = False # 正常显示负号
# Step 1: 读取附件四原始数据
def load_attachment4_data(file_path):
pref_df = pd.read_excel(file_path, sheet_name='Pref', header=None)
vwin_df = pd.read_excel(file_path, sheet_name='Vwin', header=None)
pref = pref_df.values.T # (10风机 x 300时间步)
vwin = vwin_df.values.T # (10风机 x 300时间步)
return pref, vwin
# Step 2: 卡尔曼滤波降噪
def apply_kalman(signal):
kf = KalmanFilter(transition_matrices=[1], observation_matrices=[1],
initial_state_mean=signal[0], initial_state_covariance=1,
observation_covariance=1, transition_covariance=0.01)
state_means, _ = kf.filter(signal)
return state_means.flatten()
def kalman_filter_data(data):
return np.apply_along_axis(apply_kalman, axis=1, arr=data)
# Step 3: 延迟识别与对齐
def find_lag(x, y, max_lag=10):
cross_corr = np.correlate(x - np.mean(x), y - np.mean(y), mode='full')
lags = np.arange(-max_lag, max_lag + 1)
best_lag = lags[np.argmax(cross_corr[len(cross_corr)//2 - max_lag : len(cross_corr)//2 + max_lag + 1])]
return best_lag
def align_signals(signal1, signal2, lag):
t = np.arange(len(signal1))
if lag > 0:
f = interp1d(t, signal1, kind='linear', fill_value='extrapolate')
aligned = f(t + lag)
return aligned[:len(signal2)], signal2
elif lag < 0:
f = interp1d(t, signal2, kind='linear', fill_value='extrapolate')
aligned = f(t - lag)
return signal1, aligned[:len(signal1)]
else:
return signal1, signal2
def align_all_signals(pref, vwin):
pref_aligned = np.zeros_like(pref)
vwin_aligned = np.zeros_like(vwin)
for i in range(pref.shape[0]):
lag = find_lag(pref[i], vwin[i])
p, v = align_signals(pref[i], vwin[i], lag)
pref_aligned[i] = p
vwin_aligned[i] = v
return pref_aligned, vwin_aligned
# Step 4: 计算主轴扭矩和塔架推力
def calculate_shaft_torque(power, wind_speed, rotor_radius=50, efficiency=0.9):
tip_speed_ratio = 8
angular_velocity = tip_speed_ratio * wind_speed / rotor_radius
torque = power / angular_velocity
return torque * efficiency
def calculate_tower_thrust(wind_speed, rotor_radius=50):
swept_area = np.pi * rotor_radius ** 2
thrust = 0.5 * 1.225 * swept_area * wind_speed ** 2 * 0.7
return thrust
# Step 5: 遗传算法优化功率分配(带疲劳损伤最小化 + 误差惩罚)
def optimize_power_with_ga(pref_aligned, vwin_aligned):
total_power = np.sum(pref_aligned, axis=0)
optimized_power = np.zeros_like(pref_aligned)
creator.create("FitnessMin", base.Fitness, weights=(-1.0,))
creator.create("Individual", list, fitness=creator.FitnessMin)
toolbox = base.Toolbox()
bounds = [(0, 5e6) for _ in range(10)] # 每台风机功率上限为5MW
toolbox.register("power_gene", np.random.uniform, 0, 1)
toolbox.register("individual", tools.initRepeat, creator.Individual, toolbox.power_gene, n=10)
toolbox.register("population", tools.initRepeat, list, toolbox.individual)
def eval_func(individual, t, prev_power=None):
powers = np.clip(individual, 0, 5e6)
if abs(np.sum(powers) - total_power[t]) > 1e4:
return 1e6,
wind_speeds = vwin_aligned[:, t]
# 计算载荷
torques = [calculate_shaft_torque(p, wind_speeds[i]) for i, p in enumerate(powers)]
thrusts = [calculate_tower_thrust(w) for w in wind_speeds]
# 疲劳损伤
shaft_damage = np.mean([(abs(t) / 1e8)**3 for t in torques])
tower_damage = np.mean([(abs(f) / 1.2e8)**4 for f in thrusts])
total_damage = 0.6 * shaft_damage + 0.4 * tower_damage
# 功率误差惩罚
power_error = np.abs(powers - pref_aligned[:, t])
penalty = 1e3 * np.sum(power_error > 0.1 * np.abs(pref_aligned[:, t]))
# 功率变化平滑性惩罚
smooth_penalty = 0
if prev_power is not None:
smooth_penalty = 1e2 * np.sum((powers - prev_power) ** 2)
return total_damage + penalty + smooth_penalty,
toolbox.register("evaluate", lambda ind, t, prev=None: eval_func(ind, t, prev))
toolbox.register("mate", tools.cxSimulatedBinaryBounded,
low=[b[0] for b in bounds], up=[b[1] for b in bounds], eta=20.0)
toolbox.register("mutate", tools.mutPolynomialBounded,
low=[b[0] for b in bounds], up=[b[1] for b in bounds], eta=20.0, indpb=0.1)
toolbox.register("select", tools.selTournament, tournsize=3)
prev_power = None
for t in range(300):
pop = toolbox.population(n=50)
for gen in range(10):
offspring = algorithms.varAnd(pop, toolbox, cxpb=0.7, mutpb=0.3)
fits = [toolbox.evaluate(ind, t, prev_power) for ind in offspring]
fits = [f[0] for f in fits]
selected = tools.selBest(pop + offspring, k=50)
pop = selected
best = tools.selBest(pop, k=1)[0]
optimized_power[:, t] = np.clip(best, 0, 5e6)
prev_power = optimized_power[:, t]
return optimized_power
# Step 6: 计算疲劳损伤
def compute_cumulative_damage(power, wind_speed):
cumulative_damage = []
for t in range(power.shape[1]):
torques = [calculate_shaft_torque(power[i, t], wind_speed[i, t]) for i in range(10)]
thrusts = [calculate_tower_thrust(wind_speed[i, t]) for i in range(10)]
shaft_damage = np.mean([(abs(t) / 1e8)**3 for t in torques])
tower_damage = np.mean([(abs(f) / 1.2e8)**4 for f in thrusts])
total = 0.6 * shaft_damage + 0.4 * tower_damage
cumulative_damage.append(total if t == 0 else cumulative_damage[-1] + total)
return np.array(cumulative_damage)
# Step 7: 动图可视化
def animate_damage_comparison(cumulative_damage_optimized, cumulative_damage_original):
fig, ax = plt.subplots()
line1, = ax.plot([], [], label='优化后累计损伤')
line2, = ax.plot([], [], label='原始累计损伤')
ax.set_xlim(0, 300)
ax.set_ylim(0, 2)
ax.set_xlabel('时间步')
ax.set_ylabel('累计疲劳损伤')
ax.legend()
def update(frame):
line1.set_data(range(frame), cumulative_damage_optimized[:frame])
line2.set_data(range(frame), cumulative_damage_original[:frame])
return line1, line2
ani = FuncAnimation(fig, update, frames=300, interval=100, blit=True)
plt.title('优化前后累计疲劳损伤对比')
plt.show()
# Step 8: 对比图
def plot_damage_comparison(cumulative_damage_optimized, cumulative_damage_original):
plt.plot(cumulative_damage_optimized, label='优化后累计损伤')
plt.plot(cumulative_damage_original, label='原始累计损伤')
plt.xlabel('时间步')
plt.ylabel('累计疲劳损伤')
plt.title('有无优化器的累计疲劳损伤对比')
plt.legend()
plt.grid()
plt.show()
# 主程序入口
if __name__ == "__main__":
file_path = r'C:\Users\1\Desktop\附件4-噪声和延迟作用下的采集数据.xlsx'
pref, vwin = load_attachment4_data(file_path)
# Step 2: 卡尔曼滤波降噪
pref_clean = kalman_filter_data(pref)
vwin_clean = kalman_filter_data(vwin)
# Step 3: 延迟对齐
pref_aligned, vwin_aligned = align_all_signals(pref_clean, vwin_clean)
# Step 4: 使用问题四更新模型优化功率分配
optimized_power = optimize_power_with_ga(pref_aligned, vwin_aligned)
# Step 6: 计算疲劳损伤
cumulative_damage_optimized = compute_cumulative_damage(optimized_power, vwin_aligned)
cumulative_damage_original = compute_cumulative_damage(pref_aligned, vwin_aligned)
# Step 7 & 8: 可视化
animate_damage_comparison(cumulative_damage_optimized, cumulative_damage_original)
plot_damage_comparison(cumulative_damage_optimized, cumulative_damage_original)