你希望在当前的 **遗传算法优化模型中加入新的约束条件**，即：
- 对于任意时间步 $ t $ 和风机 $ i $：
- $ |v_{i,t}^{\text{measured}} - v_{i,t}^{\text{true}}| \leq 0.1 \cdot v_{i,t}^{\text{true}} $
- $ |p_{i,t}^{\text{measured}} - p_{i,t}^{\text{true}}| \leq 0.1 \cdot p_{i,t}^{\text{true}} $
这表示**功率和风速的测量值与真实值之间的误差不超过 ±10%**。
---
## ✅ 解决思路
### 1. **新增数据输入：真实功率与真实风速**
你需要读取两个新数据源：
- `true_power_data.xlsx`：真实功率值（用于判断误差）
- `true_wind_speed_data.xlsx`：真实风速值（用于判断误差）
### 2. **在遗传算法的目标函数中添加误差约束惩罚项**
如果某个风机的测量值与真实值的误差超过 ±10%，则在目标函数中引入惩罚项，使其适应度值变差。
---
## ✅ 修改后的完整代码（新增约束）
```python
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
# ================== 参数设置 ==================
input_damage_file = '附件6-问题二答案表.xlsx' # 来自问题二的输出文件
input_power_file = 'output_data.xlsx' # 原始功率数据文件
input_true_power_file = 'true_power_data.xlsx' # 新增：真实功率数据
input_true_wind_file = 'true_wind_speed_data.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
# ================== 读取问题二输出的载荷数据 ==================
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 # (100风机 x 2000时间步)
tower_df = pd.read_excel(file_path, sheet_name='塔架推力', header=0)
tower_forces = tower_df.iloc[:, 1:].values.T # (100风机 x 2000时间步)
return shaft_torques, tower_forces
# ================== 读取原始功率数据 ==================
def load_power_data(file_path, sheet_name='WF_1'):
df = pd.read_excel(file_path, sheet_name=sheet_name, header=None)
data = df.values
num_turbines = 100
num_time_steps = 2000
powers = []
for i in range(num_turbines):
start = i * num_time_steps
end = start + num_time_steps
powers.append(data[start:end, 1]) # 第2列是功率
return np.array(powers) # (100风机 x 2000时间步)
# ================== 新增：读取真实功率与风速数据 ==================
def load_true_data(power_path, wind_path):
true_power = pd.read_excel(power_path, header=None).values
true_wind = pd.read_excel(wind_path, header=None).values
num_turbines = 100
num_time_steps = 2000
true_power_list = []
true_wind_list = []
for i in range(num_turbines):
start = i * num_time_steps
end = start + num_time_steps
true_power_list.append(true_power[start:end, 1])
true_wind_list.append(true_wind[start:end, 1])
return np.array(true_power_list), np.array(true_wind_list) # (100 x 2000)
# ================== 疲劳损伤计算函数 ==================
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):
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.mean(np.abs(powers - true_power) / (np.abs(true_power) + 1e-6))
wind_error = np.mean(np.abs(powers - true_wind) / (np.abs(true_wind) + 1e-6))
# 如果误差超过10%，引入惩罚
power_penalty = 1e3 * np.sum((np.abs(powers - true_power) > 0.1 * np.abs(true_power)))
wind_penalty = 1e3 * np.sum((np.abs(powers - true_wind) > 0.1 * np.abs(true_wind)))
# 总目标函数 = 损伤 + 惩罚项
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):
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)
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)[0]
return best_powers, best_damage
# ================== 主程序入口 ==================
if __name__ == "__main__":
# Step 1: 加载问题二的数据
shaft_torques, tower_forces = load_damage_data(input_damage_file)
power_data = load_power_data(input_power_file)
# Step 2: 加载真实功率与风速数据
true_power_data, true_wind_data = load_true_data(input_true_power_file, input_true_wind_file)
n_time_steps = 2000
optimized_results = []
# 存储历史数据用于动画展示
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()
ani = None # 定义全局动画对象
def run_animation(t):
print(f"处理第 {t + 1} 时间步...")
shaft_t = shaft_torques[:, t]
tower_f = tower_forces[:, t]
power_initial = power_data[:, t]
# 获取真实值
true_power = true_power_data[:, t]
true_wind = true_wind_data[:, t]
opt_powers, damage_after = optimize_ga(shaft_t, tower_f, power_initial, true_power, true_wind)
total_power_after = np.sum(opt_powers)
constraint_ok = abs(np.sum(opt_powers) - np.sum(power_initial)) < 1e4
# 计算损伤（优化前 vs 优化后）
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_power, true_wind)[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(n_time_steps):
run_animation(t)
# 保存动画为 GIF
ani = animation.FuncAnimation(fig, animate, frames=n_time_steps, interval=1000, repeat=False)
ani.save('real_time_optimization_with_noise_constraint.gif', writer='pillow')
plt.close()
# 将结果转换为 numpy 数组 (2000时间步 x 100风机)
optimized_array = np.array(optimized_results) # shape: (2000, 100)
# 添加时间列
n_time_steps = 2000
time_column = np.arange(1, n_time_steps + 1).reshape(-1, 1) # shape: (2000, 1)
output_data = np.hstack([time_column, optimized_array]) # shape: (2000, 101)
# 构建 DataFrame
columns = ['时间'] + [f'风机{i+1}' for i in range(100)]
result_df = pd.DataFrame(output_data, columns=columns)
# 写入 Excel 文件
result_df.to_excel(output_result_file, index=False)
print(f"优化完成，结果已保存至 {output_result_file}")
```
---
## ✅ 新增功能说明
| 模块 | 功能 |
|------|------|
| `load_true_data` | 读取真实功率和风速数据 |
| `objective_function` | 新增误差惩罚项，确保测量值与真实值误差在 ±10% 内 |
| `optimize_ga` | 传入真实功率与风速参数 |
| `动画显示` | 保留实时动画展示优化过程 |
---
###
最低0.47元/天 解锁文章
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
1553
