2024年中国研究生数学建模竞赛A题 风电场有功功率优化分配(空气动力学机理建模+鲁棒优化)

2024年中国研究生数学建模竞赛A题 风电场有功功率优化分配(空气动力学机理建模+鲁棒优化)

已完成更新，代码量近千行，进阶可视化近10张，篇幅有限，此处仅放出部分内容~

问题分析

问题背景和重述：

随着我国风电行业的迅速发展，大型风力发电机组和大规模风电场逐步投入使用。然而，大型风机的机械部件更加柔性，导致其疲劳损伤累积速度加快。这不仅增加了风机的维护成本，降低了发电效率，还可能在极端情况下引发风机倒塌、叶片断裂等安全风险。

因此，亟需通过优化手段，降低风机运行过程中的累积疲劳损伤，减少因疲劳导致的寿命缩短风险。降低风机的疲劳损伤能够显著减少机械系统的故障率，为风电场运营商带来巨大的经济和社会效益。

问题描述：

在风电场的运行过程中，电网会下发有功功率调度指令 $P_t$ 到场站控制端。场站的自动发电控制系统（AGC）以每秒钟一次的频率响应调度指令，计算场内每台风机需要发出的功率，并将结果发送给各风机。传统的平均分配方法将 $P_t$ 平均分配给每台风机，但未考虑风机的运行状态和累积疲劳损伤之间的关系，导致部分风机长期处于高疲劳负荷状态。

需要解决的问题：

1. 问题一： 建立一个低复杂度的数学模型，实时计算风机主轴和塔架的累积疲劳损伤程度。要求计算速度快（小于1秒），结果能正确反映疲劳损伤程度，并与雨流计数法的结果具有相似性。

2. 问题二： 建立数学模型，根据风速和功率参考值，估算风机所承受的塔架推力和主轴扭矩。要求模型合理有效，能够准确估算应力/扭矩值，并与参考数据进行对比。

3. 问题三： 建立有功功率的优化分配模型，以降低风电场所有风机的总体疲劳损伤程度。需要满足功率平衡和各风机功率限幅等约束条件，且优化计算需要每秒进行一次。

4. 问题四： 考虑通信延迟和测量噪声对优化模型的影响，完善问题三中的优化方案。要求模型对噪声和延迟具有鲁棒性，并展示优化效果的对比。

整体问题分析：

1. 风机疲劳损伤的计算挑战：

实时性要求高： 需要在每秒内完成计算，以满足AGC系统的实时控制需求。
数据复杂性： 载荷数据随机性强，传统的雨流计数法计算复杂，无法实时应用。
模型简化需求： 需要在保证计算精度的前提下，简化模型以降低计算复杂度。

2. 风机受力估算的难点：

物理模型复杂： 风机的受力情况与风速、功率输出等因素密切相关，且关系非线性。
数据有限： 需要在有限的数据条件下，建立准确的受力估算模型。

3. 优化模型的构建与求解挑战：

多目标优化： 需要同时考虑主轴和塔架的疲劳损伤。
高维度： 风电场包含大量风机，优化变量数量大。
实时性： 优化计算需要在每秒内完成，要求算法高效。

4. 噪声和延迟的影响：

测量噪声： 传感器数据存在随机误差，影响模型输入的准确性。
通信延迟： 数据传输过程中存在延迟，导致部分数据无法及时获取。

数据说明分析：

附件1（疲劳评估数据）：

包含100台风机在100秒内的主轴扭矩和塔架推力数据。
提供了基于雨流计数法计算的累积疲劳损伤值和等效疲劳载荷，可用于验证模型的准确性。

附件2（风电机组采集数据）：

包含两个风电场，每个风电场有100台风机的数据。
数据包括时间、输入（功率指令、风速）、输出（扭矩、推力、实际功率）和状态（桨距角、转速）等。
可用于建立风机受力估算模型和优化模型。

附件3（噪声和延迟作用下的采集数据）：

在附件2的基础上添加了测量噪声和通信延迟，模拟了实际运行环境。
用于测试优化模型在噪声和延迟条件下的鲁棒性。

题目思路

问题一的求解思路：

目标： 建立一个低复杂度的模型，实时计算风机主轴和塔架的累积疲劳损伤。

思路：

1. 简化疲劳损伤计算：

   采用应力幅值统计法： 使用应力的统计特征（如均方根值、标准差、峰值）来近似评估疲劳损伤。
   等效应力法： 将复杂的随机应力序列转换为等效的恒幅应力。

2. 基于S-N曲线的计算：

   使用材料的S-N曲线，将等效应力映射到疲劳寿命（循环次数 $N_f$）。
   计算每秒的疲劳损伤增量 $D=\frac{\Delta n}{N_f}$，其中 $\Delta n=1$（每秒一个循环假设）。

3. 累积疲劳损伤：

   对每秒的疲劳损伤增量进行累加，得到累积疲劳损伤 $D_{\text{total}}=\sum D$。

4. 验证模型：

   将计算结果与附件1中雨流计数法的结果进行比较，评估模型的准确性。
   选取5-10个代表性的风机，绘制累积疲劳损伤随时间的增长曲线。

问题二的求解思路：

目标： 根据风速和功率参考值，估算风机的塔架推力和主轴扭矩。

思路：

1. 建立物理模型：

   风轮动力学模型： 基于空气动力学原理，利用风轮的功率和受力公式。
   功率系数 $C_p$ 和推力系数 $C_t$： 反映风机在不同工况下的性能。

2. 公式推导：

   功率公式：
   $ P = \frac{1}{2} \rho A C_p V^3 $
   其中 $\rho$ 为空气密度，$A$ 为扫风面积，$V$ 为风速。

   推力公式：
   $$F_t=\frac{1}{2}\rho A{C}_t{V}^2$$

   扭矩公式：
   $$T_s=\frac{P}{\omega }$$
   其中 $\omega$ 为转速。

3. 参数获取：

   系数拟合： 利用附件2的数据，拟合出 $C_p$ 和 $C_t$ 随桨距角、风速等参数的变化关系。

4. 模型计算：

   输入风速 $V$ 和功率参考值 $P_{\text{ref}}$，计算对应的 $F_t$ 和$T_s$。

5. 误差分析：

   计算估算值与参考值的误差平方和，评估模型的准确性。

问题三的求解思路：

目标： 建立有功功率优化分配模型，降低整体疲劳损伤。

思路：

1. 定义优化目标函数：

   多目标优化： 最小化所有风机主轴和塔架的累积疲劳损伤。
   $$\min \left(\sum _{i=1}^N{D}_{\text{shaft}}^i+\sum _{i=1}^N{D}_{\text{tower}}^i\right)$$

2. 约束条件：

   功率平衡：
   $$\sum _{i=1}^N{P}_{\text{ref}}^i=P_t$$

   功率限幅：
   $$0\le P_{\text{ref}}^i\le P_{\text{rated}}$$

   偏差限制：
   $$|P_{\text{ref}}^i{\bar{P}}_{\text{ref}}^i|\le 1\text{MW}$$
   其中 $${\bar{P}}_{\text{ref}}^i=\frac{P_t}{N}$$。

3. 优化变量：

   风机的功率参考值 $P_{\text{ref}}^i$，$i=1,2,...,N$。

4. 优化算法：

   线性或非线性规划： 根据目标函数和约束的性质，选择合适的优化算法。
   算法效率： 需要保证算法在1秒内完成计算。

5. 模型验证：

   对比优化前后的疲劳损伤累积值、功率分配的均衡性。
   绘制实时计算结果的动图，展示优化过程。

问题四的求解思路：

目标： 考虑测量噪声和通信延迟，提升优化模型的鲁棒性。

思路：

1. 处理测量噪声：

   数据滤波： 使用卡尔曼滤波或移动平均滤波，对测量数据进行平滑。
   鲁棒优化： 在模型中考虑噪声的影响，采用鲁棒优化方法。

2. 处理通信延迟：

   数据预测： 利用历史数据，预测当前的风速和状态参数。
   延迟补偿： 在优化模型中引入延迟项，调整优化策略。

3. 模型调整：

   调整目标函数： 考虑噪声和延迟对疲劳损伤的影响，重新定义目标函数。
   调整约束条件： 放宽或修改部分约束，以适应数据的不确定性。

4. 模型验证：

   使用附件4的数据，测试模型在噪声和延迟条件下的性能。
   比较有无鲁棒性处理时的优化效果，展示模型对噪声和延迟的抑制能力。

总结：

通过上述分析，我们需要建立一个能够实时计算风机累积疲劳损伤的简化模型，准确估算风机的受力情况，并在此基础上构建一个高效的有功功率优化分配模型。模型需要考虑实际运行中的测量噪声和通信延迟，具备良好的鲁棒性。最终的目标是降低风电场的整体疲劳损伤，延长风机的使用寿命，提高风电场的经济效益和安全性。

求解过程

问题一：

引入了等效疲劳载荷的概念，使用实时的递推算法，降低了计算复杂度。
通过材料的 S-N 曲线，直接计算累积疲劳损伤增量，避免了复杂的循环计数。

问题二：

采用数据驱动的方法，利用实际数据拟合 $C_p$ 和 $C_t$ 的关系，提高了模型的准确性。
简化了空气动力学模型，使其适用于实时计算。

问题三：

引入了自适应权重的优化目标函数，根据疲劳损伤的实时变化调整优化方向。
采用高效的优化算法，保证了实时性的要求。

问题四：

利用卡尔曼滤波器和鲁棒优化方法，增强了模型对噪声和延迟的抵抗能力。
创新性地将不确定性引入优化模型，使得优化结果更加可靠。

问题一：风机主轴及塔架疲劳损伤程度量化指标计算低复杂度模型

建模过程：

为了实时计算风机主轴和塔架的累积疲劳损伤，我们需要一个计算速度快、精度较高的模型。传统的雨流计数法虽然精度高，但计算复杂度过高，无法满足实时性要求。为此，我们可以采用基于应力幅值的简化模型。

1. 等效疲劳载荷计算

等效疲劳载荷（Equivalent Fatigue Load，EFL）是将随机载荷序列转换为等效的恒幅载荷，以便于疲劳损伤计算。EFL可以通过以下公式计算：

$$EFL=(\frac{1}{T}\sum _{i=1}^n{F}_i^m\Delta t_i{)}^{1/m}$$

其中：

$F_i$ 为载荷在第 $i$ 个时间间隔的瞬时值。
$\Delta t_i$ 为时间间隔，取1秒。
$T$ 为总时间。
$m$ 为材料的 Wohler 曲线斜率（S-N 曲线中的疲劳强度指数）。

创新点：

实时性： 通过引入滑动窗口或递推公式，避免对全局数据的重复计算。
低复杂度： 使用简化的统计方法，如均方根值（RMS）等。

2. 累积疲劳损伤计算

根据 Palmgren-Miner 线性累积损伤理论，累积疲劳损伤可以表示为：

$$D=\frac{n}{N_f}$$

其中：

$n$ 为实际循环次数。
$N_f$ 为在等效应力幅下的疲劳寿命。

由于实时计算，每秒可以视为一个循环，即 $n=1$。那么累积损伤增量为：

$$省略部分内容$$

3. 计算流程：

步骤1： 采集实时载荷数据 $F(t)$。
步骤2： 计算瞬时等效载荷 $EFL$。
步骤3： 利用材料的 S-N 曲线，查找对应的疲劳寿命$N_f$。
步骤4： 计算累积损伤增量 $\Delta D=\frac{1}{N_f}$。
步骤5： 更新累积疲劳损伤 $D_{\text{total}}=D_{\text{total}}+\Delta D$。

示例代码：
![在这里插入图片描述](https://i-blog.csdnimg.cn/direct/77f048e82189402e9a1a28db511f2b63.png

import numpy as np
import pandas as pd

# 假设 m = 3 （材料参数），T = 1s
m = 3

# 读取载荷数据（以主轴扭矩为例）
torque_data = pd.read_excel('附件1-疲劳评估数据.xls', sheet_name='主轴扭矩')

# 初始化累积损伤
D_total = np.zeros(len(torque_data.columns) 1)

# 遍历时间序列
for t in range(1, len(torque_data)):
    # 获取当前时刻所有风机的载荷
    F_t = torque_data.iloc[t, 1:].values
    # 计算等效载荷
    EFL = F_t ** m
    # 假设 S-N 曲线为 N_f = K / EFL^m，K 为材料常数
    K = 1e12  # 假设值，需要根据实际材料参数确定
    N_f = K / EFL

    # 省略部分内容


# 输出结果
print('累积疲劳损伤：', D_total)

模型验证：

将计算结果与附件1中提供的累积疲劳损伤值进行对比，选取5-10个风机，绘制累积疲劳损伤随时间的增长曲线，观察趋势是否一致。

问题二：利用风速及功率估算塔架推力和主轴扭矩

建模过程：

为了估算塔架推力和主轴扭矩，我们需要建立风机的空气动力学模型。

1. 风机的空气动力学模型

功率方程：

$$P=\frac{1}{2}\rho A{C}_p(\lambda ,\beta )V^3$$

推力方程：

$$F_t=\frac{1}{2}\rho A{C}_t(\lambda ,\beta )V^2$$

其中：

$\rho$ 为空气密度。
$A$ 为风轮扫风面积。
$C_p$ 为功率系数。
$C_t$ 为推力系数。
$\lambda$ 为速比，$\lambda =\frac{\omega R}{V}$。
$\beta$ 为桨距角。
$V$ 为风速。
$R$ 为风轮半径。
$\omega$ 为角速度。

创新点：

简化模型： 将 $C_p$ 和 $C_t$ 表示为速比和桨距角的函数，可以采用拟合或查表的方法。
数据驱动建模： 利用附件2的数据，拟合 $C_p$ 和 $C_t$ 的表达式，提高模型精度。

2. 模型建立

步骤1： 计算速比$\lambda$。

$$\lambda =\frac{\omega R}{V}$$

步骤2： 拟合 $C_p$ 和 $C_t$。

使用多元线性回归或机器学习方法（如随机森林、神经网络）拟合 $C_p(\lambda ,\beta )$ 和 $C_t(\lambda ,\beta )$。

步骤3： 估算推力和扭矩。

$$T_s=\frac{P}{\omega }$$

$$F_t=\frac{1}{2}\rho A{C}_t{V}^2$$

3. 计算流程：

输入： 风速 $V$、功率参考值 $P_{\text{ref}}$、桨距角 $\beta$、转速 $\omega$。
步骤1： 计算速比 $\lambda$。
步骤2： 计算 $C_p$ 和 $C_t$。
步骤3： 计算推力$F_t$ 和扭矩 $T_s$。

示例代码：
![请添加图片描述](https://i-blog.csdnimg.cn/direct/f834376ffb8e4c1ab29f158c8b6cea6d.png

import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression

# 读取数据
data = pd.read_excel('附件2-风电机组采集数据.xlsx')

# 拟合 C_p 和 C_t
def fit_coefficients(data):
    # 提取特征和目标变量
    X = data[['lambda', 'beta']]
    y_cp = data['Cp']
    y_ct = data['Ct']
    # 拟合模型

    # 省略部分内容

    return model_cp, model_ct

# 计算速比
def calculate_lambda(omega, R, V):
    return omega * R / V

# 主函数
def estimate_Ts_Ft(V, P_ref, beta, omega, R, rho, A):
    # 计算速比
    lambd = calculate_lambda(omega, R, V)
    # 拟合 C_p 和 C_t
    model_cp, model_ct = fit_coefficients(data)
    # 预测 C_p 和 C_t

    # 省略部分内容

    # 省略部分内容

    return Ts, Ft

4. 模型验证：

计算所有时刻的 $T_s$ 和 $F_t$，与附件2中给出的参考值进行对比，计算误差平方和。

问题三：有功调度优化问题构建与实时求解

建模过程：

需要构建一个实时的有功功率优化分配模型，以降低风电场所有风机的累积疲劳损伤。

1. 优化目标函数

目标1： 最小化主轴累积疲劳损伤总和 $D_{\text{shaft}}$。
目标2： 最小化塔架累积疲劳损伤总和 $D_{\text{tower}}$。

可以将两个目标合并为一个加权和：

$$\min \left(w_1\sum _{i=1}^N{D}_{\text{shaft}}^i+w_2\sum _{i=1}^N{D}_{\text{tower}}^i\right)$$

其中 $w_1$ 和 $w_2$ 为权重，可以根据实际需求设定。

创新点：

引入权重自适应调整： 根据实时的疲劳损伤程度，动态调整权重 $w_1$ 和 $w_2$，使得优化更加智能化。
采用预测模型： 结合问题二的模型，预测在不同功率分配下的疲劳损伤。

2. 约束条件

功率平衡：

$$\sum _{i=1}^N{P}_{\text{ref}}^i=P_t$$

功率限幅：

$$0\le P_{\text{ref}}^i\le P_{\text{rated}}$$

偏差限制：

$$|P_{\text{ref}}^i{\bar{P}}_{\text{ref}}^i|\le 1\text{MW}$$

其中 ${\bar{P}}_{\text{ref}}^i=\frac{P_t}{N}$。

3. 优化变量

风机的功率参考值 $P_{\text{ref}}^i$。

4. 优化求解

算法选择：

由于目标函数和约束条件都是线性的或可以线性化，采用线性规划（Linear Programming，LP）求解。
如果目标函数是非线性的，可以采用二次规划（Quadratic Programming，QP）或其他非线性优化算法。

实时性保障：

省略部分内容

5. 计算流程：

步骤1： 收集当前风速 $V_w^i$ 和初始疲劳损伤 $D_{\text{shaft}}^i$、$D_{\text{tower}}^i$。
步骤2： 利用问题二的模型，预测在不同 $P_{\text{ref}}^i$ 下的疲劳损伤增量。
步骤3： 建立优化模型，定义目标函数和约束条件。
步骤4： 求解优化问题，得到最优的 $P_{\text{ref}}^i$。
步骤5： 将结果发送至各风机。

示例代码：

import numpy as np
from scipy.optimize import minimize

# 风机数量
N = 100

# 输入数据
P_t = ...  # 总调度指令
P_rated = 5  # 额定功率

# 初始疲劳损伤（假设已计算）
D_shaft = np.zeros(N)
D_tower = np.zeros(N)

# 平均功率分配
P_ref_avg = P_t / N

# 目标函数
def objective(P_ref):
    # 省略部分内容

# 约束条件
constraints = [
    {'type': 'eq', 'fun': lambda P_ref: np.sum(P_ref) P_t},
    {'type': 'ineq', 'fun': lambda P_ref: P_rated P_ref},
    {'type': 'ineq', 'fun': lambda P_ref: P_ref},
    {'type': 'ineq', 'fun': lambda P_ref: 1 np.abs(P_ref P_ref_avg)}
]

# 初始猜测
P_ref0 = np.full(N, P_ref_avg)

# 求解优化问题
result = minimize(objective, P_ref0, constraints=constraints, method='SLSQP')

# 最优功率分配
P_ref_opt = result.x

6. 结果分析：

对比优化前后的累积疲劳损伤。
计算功率分配的方差，评估功率分配的均衡性。

问题四：考虑通信延迟和测量噪声的有功功率优化与求解

建模过程：

在实际应用中，需要考虑测量噪声和通信延迟对优化模型的影响。

1. 处理测量噪声

数据滤波：

采用卡尔曼滤波器，对测量数据进行估计和滤波。
建立状态空间模型，估计真实的风速和疲劳损伤。

创新点：

自适应滤波器： 根据噪声水平，动态调整滤波器参数，提高滤波效果。

2. 处理通信延迟

数据预测：

利用历史数据，采用时间序列预测方法（如 ARIMA、LSTM）预测当前的风速和状态。

延迟补偿：

在优化模型中引入延迟项，使用上一次有效数据进行优化。

创新点：

鲁棒优化：

将测量噪声和延迟视为不确定性，采用鲁棒优化方法，保证在最坏情况下的优化效果。

3. 调整优化模型

新的目标函数：

考虑噪声和延迟的影响，加入对不确定性的惩罚项。

$$\min \left(w_1\sum _{i=1}^N{\hat{D}}_{\text{shaft}}^i+w_2\sum _{i=1}^N{\hat{D}}_{\text{tower}}^i+w_3\cdot \text{UncertaintyPenalty}\right)$$

约束条件：

考虑延迟情况下的约束，允许一定的功率偏差。

4. 优化求解

采用鲁棒优化算法：

使用场景规划或随机规划方法，考虑不确定性。

实时性保障：

由于鲁棒优化计算复杂度较高，可以采用近似算法或分布式计算。

示例代码：

import numpy as np
from scipy.optimize import minimize

# 假设测量噪声为高斯分布，通信延迟为 tau 秒

# 滤波器初始化
def kalman_filter(z, x_prev, P_prev, Q, R):
    # 预测
    x_pred = x_prev
    P_pred = P_prev + Q
    # 更新
    K = P_pred / (P_pred + R)
    x_updated = x_pred + K * (z x_pred)
    P_updated = (1 K) * P_pred
    return x_updated, P_updated

# 调整后的目标函数
def robust_objective(P_ref):
    # 估计的疲劳损伤增量，考虑不确定性
    delta_D_shaft = ...  # 需要考虑噪声和延迟
    delta_D_tower = ...
    # 不确定性惩罚项
    uncertainty_penalty = ...  # 根据估计误差计算
    # 总目标
    total_damage = w1 * np.sum(D_shaft + delta_D_shaft) + w2 * np.sum(D_tower + delta_D_tower) + w3 * uncertainty_penalty
    return total_damage

# 约束条件与之前类似，但需要考虑延迟

# 求解优化问题
result = minimize(robust_objective, P_ref0, constraints=constraints, method='SLSQP')

5. 模型验证：

使用附件4的数据，测试模型在噪声和延迟条件下的性能。
对比有无鲁棒优化情况下的结果，展示模型对噪声和延迟的抑制能力。

整体建模流程的代码集成示例：

通过深入挖掘和细化建模过程，我们建立了一个兼顾实时性和准确性的风电场有功功率优化分配模型。该模型创新性地结合了疲劳损伤的实时估计、风机空气动力学的简化建模、自适应的优化目标函数，以及对测量噪声和通信延迟的鲁棒处理。通过示例代码的实现，为实际工程应用提供了可行的解决方案。