import numpy as np
import pandas as pd
from deap import base, tools, algorithms, creator
import matplotlib.pyplot as plt
import matplotlib.animation as animation
# 设置中文字体支持
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
# ================== 参数设置 ==================
input_damage_file = 'C:/Users/1/Desktop/附件6-问题二答案表.xlsx'
output_result_file = '优化后的功率分配_问题三加权损伤_含噪声约束.xlsx'
# S-N曲线参数
S0_shaft = 1e8
m_shaft = 3
S0_tower = 1.2e8
m_tower = 4
# 权重
weight_shaft = 0.6
weight_tower = 0.4
# 新增文件路径
true_file = r'C:/Users/1/Desktop/问题四_延时补偿后数据.xlsx'
measured_file = r'C:/Users/1/Desktop/附件3-噪声和延迟作用下的采集数据.xlsx'
# 每台风机数据行数
rows_per_turbine = 2000
num_turbines = 100
num_time_steps = 2000
# ================== 读取问题二输出的载荷数据 ==================
def load_damage_data(file_path):
shaft_df = pd.read_excel(file_path, sheet_name='主轴扭矩', header=0)
shaft_torques = shaft_df.iloc[:, 1:].values.T
tower_df = pd.read_excel(file_path, sheet_name='塔架推力', header=0)
tower_forces = tower_df.iloc[:, 1:].values.T
return shaft_torques, tower_forces
# ================== 新增：读取真实功率与风速 + 测量功率与风速 ==================
def load_true_and_measured_data(true_path, measured_path):
true_data = pd.read_excel(true_path, header=None).values
measured_data = pd.read_excel(measured_path, header=None).values
true_power_list = []
true_wind_list = []
measured_power_list = []
measured_wind_list = []
for i in range(num_turbines):
start = i * rows_per_turbine
end = start + rows_per_turbine
true_power_list.append(true_data[start:end, 1])
true_wind_list.append(true_data[start:end, 2])
measured_power_list.append(measured_data[start:end, 1])
measured_wind_list.append(measured_data[start:end, 2])
return {
'true_power': np.array(true_power_list),
'true_wind': np.array(true_wind_list),
'measured_power': np.array(measured_power_list),
'measured_wind': np.array(measured_wind_list)
}
# ================== 疲劳损伤计算函数 ==================
def calculate_shaft_damage(torque):
stress = abs(torque)
return (stress / S0_shaft) ** m_shaft
def calculate_tower_damage(force):
stress = abs(force)
return (stress / S0_tower) ** m_tower
# ================== 目标函数（新增：噪声约束惩罚项） ==================
def objective_function(individual, shaft_torques, tower_forces, total_power,
true_power, true_wind, measured_power, measured_wind):
powers = np.clip(individual, 0, 5e6)
if abs(np.sum(powers) - total_power) > 1e4:
return 1e6,
shaft_damages = [calculate_shaft_damage(t) for t in shaft_torques]
tower_damages = [calculate_tower_damage(f) for f in tower_forces]
avg_shaft = np.mean(shaft_damages)
avg_tower = np.mean(tower_damages)
total_damage = weight_shaft * avg_shaft + weight_tower * avg_tower
power_error = np.abs(powers - measured_power) <= 0.1 * np.abs(true_power)
wind_error = np.abs(powers - measured_wind) <= 0.1 * np.abs(true_wind)
power_penalty = 1e3 * np.sum(~power_error)
wind_penalty = 1e3 * np.sum(~wind_error)
final_objective = total_damage + power_penalty + wind_penalty
return final_objective,
# ================== 功率限制函数 ==================
def get_bounds(avg_power):
low = max(0, avg_power - 1e6)
high = min(5e6, avg_power + 1e6)
return [(low, high)] * 100
# ================== 遗传算法优化器 ==================
def optimize_ga(shaft_torques, tower_forces, power_initial,
true_power, true_wind, measured_power, measured_wind):
total_power = np.sum(power_initial)
avg_power = total_power / 100
bounds = get_bounds(avg_power)
creator.create("FitnessMin", base.Fitness, weights=(-1.0,))
creator.create("Individual", list, fitness=creator.FitnessMin)
toolbox = base.Toolbox()
for i in range(100):
toolbox.register(f"attr_{i}", np.random.uniform, bounds[i][0], bounds[i][1])
toolbox.register("individual", tools.initCycle, creator.Individual,
[getattr(toolbox, f"attr_{i}") for i in range(100)], n=1)
toolbox.register("population", tools.initRepeat, list, toolbox.individual)
toolbox.register("evaluate", objective_function, shaft_torques=shaft_torques,
tower_forces=tower_forces, total_power=total_power,
true_power=true_power, true_wind=true_wind,
measured_power=measured_power, measured_wind=measured_wind)
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)
pop = toolbox.population(n=50)
hof = tools.HallOfFame(1)
stats = tools.Statistics(lambda ind: ind.fitness.values)
stats.register("avg", np.mean)
stats.register("min", np.min)
algorithms.eaSimple(pop, toolbox, cxpb=0.7, mutpb=0.3, ngen=50,
stats=stats, halloffame=hof, verbose=False)
min_bounds = np.array([b[0] for b in bounds])
max_bounds = np.array([b[1] for b in bounds])
best_powers = np.clip(hof[0], min_bounds, max_bounds)
best_damage = objective_function(best_powers, shaft_torques, tower_forces,
total_power, true_power, true_wind,
measured_power, measured_wind)[0]
return best_powers, best_damage
# ================== 主程序入口 ==================
if __name__ == "__main__":
# Step 1: 加载问题二的数据
shaft_torques, tower_forces = load_damage_data(input_damage_file)
# Step 2: 加载真实功率与风速、测量功率与风速
data_dict = load_true_and_measured_data(true_file, measured_file)
true_power = data_dict['true_power']
true_wind = data_dict['true_wind']
measured_power = data_dict['measured_power']
measured_wind = data_dict['measured_wind']
optimized_results = []
raw_power_data = [] # 用于保存原始测量功率
# 存储历史数据用于动画展示
history = {
'time_step': [],
'total_power': [],
'damage_before': [],
'damage_after': [],
'constraint_ok': []
}
print("开始逐时间步进行功率优化分配...")
fig, axs = plt.subplots(3, 1, figsize=(10, 8))
def animate(i):
if i >= len(history['time_step']):
return
axs[0].clear()
axs[1].clear()
axs[2].clear()
time_steps = history['time_step'][:i+1]
powers = history['total_power'][:i+1]
damage_before = history['damage_before'][:i+1]
damage_after = history['damage_after'][:i+1]
constraints = history['constraint_ok'][:i+1]
axs[0].plot(time_steps, powers, label='总功率', color='blue')
axs[0].set_title(f"时间步: {time_steps[-1]} | 总功率: {powers[-1]:.2f} W")
axs[0].legend()
axs[1].plot(time_steps, damage_before, label='优化前损伤', color='orange')
axs[1].plot(time_steps, damage_after, label='优化后损伤', color='green')
axs[1].set_title("主轴+塔架加权损伤")
axs[1].legend()
status = "满足" if constraints[-1] else "不满足"
color = "green" if constraints[-1] else "red"
axs[2].text(0.5, 0.5, f"约束条件: {status}", fontsize=20, ha='center', va='center', color=color)
axs[2].axis('off')
plt.tight_layout()
def run_animation(t):
print(f"处理第 {t + 1} 时间步...")
shaft_t = shaft_torques[:, t]
tower_f = tower_forces[:, t]
power_initial = measured_power[:, t]
raw_power_data.append(power_initial.copy()) # 保存原始测量功率
true_p = true_power[:, t]
true_w = true_wind[:, t]
measured_p = measured_power[:, t]
measured_w = measured_wind[:, t]
opt_powers, damage_after = optimize_ga(shaft_t, tower_f, power_initial,
true_p, true_w, measured_p, measured_w)
total_power_after = np.sum(opt_powers)
constraint_ok = abs(np.sum(opt_powers) - np.sum(power_initial)) < 1e4
damage_before_val = weight_shaft * np.mean([calculate_shaft_damage(t) for t in shaft_t]) + \
weight_tower * np.mean([calculate_tower_damage(f) for f in tower_f])
damage_after_val = objective_function(opt_powers, shaft_t, tower_f, np.sum(power_initial),
true_p, true_w, measured_p, measured_w)[0]
history['time_step'].append(t + 1)
history['total_power'].append(total_power_after)
history['damage_before'].append(damage_before_val)
history['damage_after'].append(damage_after_val)
history['constraint_ok'].append(constraint_ok)
optimized_results.append(opt_powers)
for t in range(num_time_steps):
run_animation(t)
# 保存动画为 GIF
ani = animation.FuncAnimation(fig, animate, frames=num_time_steps, interval=1000, repeat=False)
ani.save('real_time_optimization_with_noise_constraint.gif', writer='pillow')
plt.close()
# 保存结果
optimized_array = np.array(optimized_results)
raw_array = np.array(raw_power_data)
# 构建 DataFrame
time_column = np.arange(1, num_time_steps + 1).reshape(-1, 1)
output_data = np.hstack([time_column, optimized_array])
columns = ['时间'] + [f'风机{i+1}' for i in range(100)]
result_df = pd.DataFrame(output_data, columns=columns)
result_df.to_excel(output_result_file, index=False)
print(f"优化完成，结果已保存至 {output_result_file}")
# ================== 新增：绘制对比图 ==================
# 随机挑选两个风机
turbine_indices = np.random.choice(num_turbines, 2, replace=False)
# 1. 功率分配对比图
plt.figure(figsize=(14, 6))
for i, idx in enumerate(turbine_indices):
plt.subplot(2, 1, i + 1)
plt.plot(raw_array[:, idx], label='原始功率（含噪声）', alpha=0.7)
plt.plot(optimized_array[:, idx], label='优化后功率', linewidth=2)
plt.title(f'风机 {idx + 1} 功率分配对比')
plt.xlabel('时间步')
plt.ylabel('功率 (W)')
plt.legend()
plt.tight_layout()
plt.savefig('功率分配对比图.png')
plt.show()
# 2. 累计疲劳值对比图
def calculate_cumulative_damage(power_array):
damage = []
for t in range(num_time_steps):
shaft_damages = [calculate_shaft_damage(shaft_torques[i, t]) for i in range(num_turbines)]
tower_damages = [calculate_tower_damage(tower_forces[i, t]) for i in range(num_turbines)]
total = weight_shaft * np.mean(shaft_damages) + weight_tower * np.mean(tower_damages)
damage.append(total)
return np.cumsum(damage)
raw_damage = calculate_cumulative_damage(raw_array)
opt_damage = calculate_cumulative_damage(optimized_array)
plt.figure(figsize=(10, 6))
plt.plot(raw_damage, label='原始功率下的累计疲劳值', alpha=0.8)
plt.plot(opt_damage, label='优化后功率下的累计疲劳值', linewidth=2)
plt.title('降噪与延迟处理前后的累计疲劳值对比')
plt.xlabel('时间步')
plt.ylabel('累计疲劳值')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.savefig('累计疲劳值对比图.png')
plt.show()
累计疲劳的对比图中没有显示原始功率下的累计疲劳总值曲线
最新发布
Time Shaft是一个旨在提高团队效率的项目，以时间轴为核心功能，突出团队协作与沟通。项目使用GitHub进行管理，通过飞书、微信群和腾讯会议进行即时沟通。其杀手级功能是独特的时间轴模式，便于记录和查看事件。项目已进行150位测试用户的反馈收集，目前正努力达到预期的500用户目标。团队在文档、测试用例和代码标准方面做了详尽的工作，确保软件工程质量。
最低0.47元/天 解锁文章
扫一扫
1552
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
