【研赛A题成品论文】24华为杯数学建模研赛A题成品论文+可运行代码丨免费分享

2024华为杯研究生数学建模竞赛A题精品成品论文已出！

A题风电场有功功率优化分配

一、问题分析

A题是一道工程建模与优化类问题，其目的是根据题目所给的附件数据资料分析风机主轴及塔架疲劳损伤程度，以及建立优化模型求解最优有功功率分配策略。整个问题的模型主要可以分为三大块：计算疲劳损伤程度的实时模型、计算应力/扭矩值的模型、求解最优有功功率分配策略的优化模型。相对应的，也需要参赛团队有对应方面的一些知识积累和技能掌握：第一，须具备扎实的数学建模能力，能够基于风力机的物理特性建立合理的数学模型。对于问题一和二，特别需要对风能转换、力学受力分析、疲劳损伤理论（如S-N曲线、Palmgren-Miner理论）等领域有深入理解。这涉及到风速与功率、扭矩与推力、疲劳损伤的关系建模，并将这些物理关系转换为数学表达式。第二，解决问题三和四需要对多目标优化算法有深入理解，能够设计优化问题并找到最优解。同时，求解大规模实时问题需要考虑算法的计算效率和精度。第三，处理风电场的数据时，必须具备数据分析和信号处理能力，特别是对于噪声和延迟的处理。问题四中的测量噪声和通信延迟的处理需要掌握滤波技术（如卡尔曼滤波）、信号分析方法以及鲁棒控制的理论。

问题一：要求我们利用主轴扭矩和塔架推力荷载数据分别实现对风机的主轴和塔架累积疲劳损伤程度的低复杂度计算。对于题目中给出的数据，首先应使用雨流计数法统计不同幅值载荷相应的循环次数（参考[2]），然后利用附件B中的步骤，计算等效疲劳载荷和累计疲劳损伤值，这里应当对基于雨流计数法计算所得到的累积疲劳损伤值和等效疲劳载荷与给出的数据进行对比，误差应该相对较小。为了实现对雨流统计法的简化，可以考虑[8],[9]中的方法，例如参考文献[8]中提到点的选取遵循单向选取、循环进行的原则这种处理方法使得雨流计数判定只向一个方向进行避免回退取点来进行计数判定的情况的出现降低了程序实现的复杂性。此问求解较为复杂，需要计算时间小于1.00s且求解时间要尽可能地短。也可以考虑利用能量进行计算和改进。

问题二：根据风机所处位置的风速条件和功率参考值，估算当前风机所承受的应力/扭矩。结合受力分析和能量守恒，关键是理解风速、功率与主轴扭矩和塔顶推力之间的物理关系，中间参数应包括风机实际输出功率，桨距角低速轴转速，高速轴转速等。[11]涵盖了风力机的空气动力学分析，包括主轴扭矩和塔顶推力的计算方法。[10]详细介绍了风力发电的基本原理，包括风速、功率系数和风机输出功率的计算。风机的工作过程涉及风能向机械能的转化，其主轴和塔架在此过程中承受不同的力。因此，首先需要通过风速和功率的关系来计算出风机的受力情况，进而推导出主轴扭矩和塔顶推力。风机从风中提取能量，发电功率与风速的三次方成正比。通过风轮扫掠的风能转化为风机的输出功率，能量转换效率由功率系数决定，它反映了风机叶片在不同风速条件下的工作效率。通过这个关系，可以确定风速对风机输出功率的影响。主轴扭矩是通过风机发电功率和风机转速的关系计算的。主轴扭矩与发电功率和转速成反比，而转速则与风速有关。通常，风机的转速由叶尖速比（TSR）控制，它表示叶片旋转速度与风速的比例。基于此，可以通过风速和功率关系，结合叶尖速比计算主轴的扭矩。风机的塔顶推力是风对风轮施加的水平力，它与风速的平方成正比。塔顶推力直接作用于风机的塔架，风轮的受风面积、风速、空气密度以及推力系数都是影响推力的关键因素。推力系数通常根据风机设计确定，并与叶片的角度和气动特性有关。通过这些参数，可以估算出塔顶推力的大小。问题二需要通过模型计算全部时刻的应力/扭矩计算，并展示计算结果与实际数据的对比结果。

问题三：通过优化算法对风电场的有功功率进行实时分配，以降低风机的疲劳损伤并保证功率平衡。该问题可以通过构建多目标优化模型来解决，目标是最小化风机主轴和塔架的累积疲劳损伤，同时满足有功功率总和等于电网调度指令 Pt，并遵守功率参考值的约束。模型的目标函数可以设定为最小化所有风机的总体疲劳损伤程度。单个风机的总疲劳损伤包含主轴疲劳损伤和塔架疲劳损伤两部分. 模型需要满足功率平衡约束(所有风机的有功功率参考值之和应等于电网调度指令), 功率参考值约束：每台风机的功率参考值不得超过风机的额定功率（5MW）, 平滑度约束：每台风机的功率分配值与平均分配值的差异不超过1MW. 该问题可通过连续优化方法进行求解，以确保优化过程的实时性。为了加快求解速度，可以使用梯度下降法或拉格朗日乘子法进行快速迭代，确保每秒进行一次功率分配计算。

问题四：在风电场的有功功率优化调度过程中，考虑通信延迟和测量噪声对优化结果的影响。优化模型需要在不确定性条件下保持稳健性，因此在设计时必须增强模型的鲁棒性，以应对噪声和延迟带来的影响。首先是测量噪声和通信延迟的处理，可以通过构建鲁棒优化模型对噪声进行处理，采用预测控制方法（如模型预测控制，MPC[15]）来解决延迟问题。同时，噪声在优化中可以通过设置误差容忍度来建模，以确保噪声干扰不会显著影响优化结果[16]。

二、问题一模型建立与求解

2.1 问题一求解思路

风机的主轴和塔架承受着随机波动的载荷，这种载荷的循环特征不易通过简单的方法统计，通常使用雨流计数法来分析载荷的循环特性。然而，标准雨流计数法计算复杂且无法满足实时计算需求。因此，我们提出了基于三点式雨流计数法的低复杂度算法，能够在较短的时间内有效计算累积疲劳损伤。

2.2 问题一模型建立

1. 累积疲劳损伤计算

累积疲劳损伤的计算基于Palmgren-Miner线性累积损伤理论。该理论表明材料在不同应力幅值下的损伤是可累积的，其损伤程度𝐷可表示为：

其中，ni是应力幅值Si下的循环次数，Ni是该应力幅值下的疲劳寿命。累积疲劳损伤D达到1时，材料将失效。

2. S-N曲线模型

S-N曲线用于描述材料在不同应力幅值下的疲劳寿命。Wohler指数m和常数C决定了材料的疲劳特性，其表达式为：

通过双对数变形，疲劳寿命N可计算为：

在模型中，m取值10，材料常数C=9.77×10^70，通过此公式可以计算材料在特定应力幅值下的疲劳寿命。

3. Goodman修正公式

在实际应用中，由于应力循环的平均应力S_m不为零，累积疲劳损伤需进行修正。Goodman修正公式用于计算平均应力S_m作用下的等效应力幅值S_e:

其中，S_a是应力幅值，σ_b是材料的最大强度。通过此修正公式，可以得到更准确的疲劳损伤计算结果。

4. 改进的三点式雨流计数法

为了识别载荷数据中的应力循环，我们使用了改进的三点式雨流计数法。该方法通过分析载荷数据中的波峰和波谷，识别应力循环，并统计不同幅值的循环次数。

（1）波峰与波谷识别

假设载荷时间历程为一个离散时间序列X={x1,x2,...,xn} ’findpeaks 函数的目标是从序列X中识别出局部极大值(波峰)和局部极小值(波谷)。

波峰识别的核心是找到序列中满足以下条件的点x_i:

其中i为波峰所在的位置。波谷可以通过对序列取反来识别

（2）索引排序

识别出的波峰和波谷可以表示为两个集合：

为了保证极值序列的时间顺序，需要将波峰和波谷的索引I_P和I_V合并并进行升序排列，生成极值序列E={e_1,e_2,...,e_l}及其对应的索引I_E，使得：

并且极值序列E满足时间上的顺序e_1,e_2,...,e_l.

（3）三点式雨流计数法的规则

在三点式雨流计数法中，每次考察极值序列E中的三个连续点e_i,e_(i+1),e_(i+2)。通过计算它们之间的幅值变化，可以识别闭合应力循环。

幅值计算

定义幅值变化ΔS1和ΔS2为：

根据这两个幅值变化，判断是否识别到闭合循环。

闭合循环识别条件

闭合循环的条件是：

如果满足该条件，则认为e_i和e_(i+1)之间存在一个完整的应力循环，幅值为ΔS1,均值为：

此时，将e_i和e_(i+1)从极值序列中移除，并继续处理新的极值序列。

无闭合循环的处理：

如果不满足闭合循环的条件，即ΔS1>ΔS2,则移动到下一个三点组合e_(i+1),e_(i+2),e_(i+3)进行分析。

2.3 问题一模型求解与分析

为了展示模型的计算结果，选择5-10个有代表性的风机，绘制主轴扭矩和塔架推力的时序数据、等效疲劳载荷及累积疲劳损伤的累积图，并与参考数据进行对比。最后，输出所有100台风机在100秒时长内每秒的累积疲劳损伤值。

为了保证模型的实时性，我们对标准雨流计数法进行了简化，同时减少了复杂的迭代步骤。整个计算过程可以在1秒内完成，适合在CPU环境下实时运行。

综上所述，基于三点式雨流计数法、Goodman修正及S-N曲线模型的低复杂度算法能够高效、准确地计算风机主轴及塔架的累积疲劳损伤值，建模方法100台风机的所有数据样本均有效，并满足实时计算需求。
2.4 参考代码

clc, clear, close all

% ================== 读取 .xls 文件 ==================
filename = 'data1.xls';  % 文件名

% 读取第2个工作表（主轴扭矩数据）
torque_data = xlsread(filename, 2);  % 主轴扭矩
time = torque_data(1:end-2, 1);  % 第一列为时间
torque_wind_turbines = torque_data(1:end-2, 2:end);  % 每台风机的扭矩载荷时序数据

% 读取第3个工作表（塔架推力数据）
thrust_data = xlsread(filename, 3);  % 塔架推力
thrust_wind_turbines = thrust_data(1:end-2, 2:end);  % 每台风机的塔架推力载荷时序数据

% 提取参考的等效疲劳载荷与累积疲劳损伤值（最后两行）
ref_equivalent_load_torque = torque_data(end-1, 2:end);  % 主轴扭矩的等效疲劳载荷
ref_fatigue_damage_torque = torque_data(end, 2:end);  % 主轴扭矩的累积疲劳损伤值
ref_equivalent_load_thrust = thrust_data(end-1, 2:end);  % 塔架推力的等效疲劳载荷
ref_fatigue_damage_thrust = thrust_data(end, 2:end);  % 塔架推力的累积疲劳损伤值

% 材料参数
sigma_b = 50000000;  % 材料拉伸断裂时的最大载荷值
m = 10;              % Wohler指数
C = 9.77e70;         % S-N曲线常数
design_life = 42565440.4361;  % 风机设计寿命,
NPeaks = 100;

% ================== 计算每台风机的指标 ==================
num_wind_turbines = size(torque_wind_turbines, 2);  % 风机数量

calculated_equivalent_load_torque = zeros(1, num_wind_turbines);
calculated_fatigue_damage_torque = zeros(1, num_wind_turbines);
calculated_equivalent_load_thrust = zeros(1, num_wind_turbines);
calculated_fatigue_damage_thrust = zeros(1, num_wind_turbines);

cpu_times = zeros(1, num_wind_turbines);  % 记录每台风机的计算时间

for i = 1:num_wind_turbines
    % 开始计时
    tic;
    
    % 提取每台风机的主轴扭矩数据
    torque_data_turbine = torque_wind_turbines(:, i);
    
    % 提取每台风机的塔架推力数据
    thrust_data_turbine = thrust_wind_turbines(:, i);
    
    % 计算主轴扭矩的等效疲劳载荷与累积疲劳损伤值
    [ranges_torque, means_torque, counts_torque] = three_point_rainflow(torque_data_turbine, NPeaks);  % 三点式雨流计数法
    L_i_torque = arrayfun(@(r, m) goodman_correction(r, m, sigma_b), ranges_torque, means_torque);
    [calculated_equivalent_load_torque(i), equivalent_load_torque_vector] = calculate_equivalent_load(L_i_torque, counts_torque, m, design_life);
    [calculated_fatigue_damage_torque(i), fatigue_damage_torque_vector] = calculate_fatigue_damage(L_i_torque, counts_torque, C, m);
    
    % 计算塔架推力的等效疲劳载荷与累积疲劳损伤值
    [ranges_thrust, means_thrust, counts_thrust] = three_point_rainflow(thrust_data_turbine, NPeaks);  % 三点式雨流计数法
    L_i_thrust = arrayfun(@(r, m) goodman_correction(r, m, sigma_b), ranges_thrust, means_thrust);
    [calculated_equivalent_load_thrust(i), equivalent_load_thrust_vector] = calculate_equivalent_load(L_i_thrust, counts_thrust, m, design_life);
    [calculated_fatigue_damage_thrust(i), fatigue_damage_thrust_vector] = calculate_fatigue_damage(L_i_thrust, counts_thrust, C, m);
    
    % 停止计时，记录时间
    cpu_times(i) = toc;
    
    % 输出中间变量图像
    if mod(i, 25) == 0  % 每隔25个风机选取一个
        figure;
        
        subplot(4, 1, 1);
        plot(time, torque_data_turbine);
        title(sprintf('风机 %d 的主轴扭矩', i));
        xlabel('时间 (s)');
        ylabel('主轴扭矩 (Nm)');
        
        subplot(4, 1, 2);
        plot(time, thrust_data_turbine);
        title(sprintf('风机 %d 的塔架推力', i));
        xlabel('时间 (s)');
        ylabel('塔架推力 (N)');
        
        % 分开绘制等效疲劳载荷和累积疲劳损伤值的累计图
        subplot(4, 1, 3);
        plot(cumsum(equivalent_load_torque_vector), 'b', 'DisplayName', '等效疲劳载荷累计值');
        title(sprintf('风机 %d 的等效疲劳载荷累计', i));
        xlabel('索引');
        ylabel('累积等效疲劳载荷');
        legend('show');
        grid on;
        
        subplot(4, 1, 4);
        plot(cumsum(fatigue_damage_torque_vector), 'r', 'DisplayName', '累积疲劳损伤值');
        title(sprintf('风机 %d 的累积疲劳损伤值', i));
        xlabel('索引');
        ylabel('累积疲劳损伤');
        legend('show');
        grid on;
    end
end

% ================== 结果对比与图形化 ==================
% 绘制每台风机的CPU时间消耗
figure;
plot(1:num_wind_turbines, cpu_times, 'b-o', 'LineWidth', 2);
title('每台风机计算的CPU时间');
xlabel('风机编号');
ylabel('CPU时间 (秒)');
grid on;

% 主轴扭矩的等效疲劳载荷与参考值对比
figure;
subplot(2,1,1);
plot(1:num_wind_turbines, calculated_equivalent_load_torque, 'b', 'DisplayName', '计算值');
hold on;
plot(1:num_wind_turbines, ref_equivalent_load_torque, 'r--', 'DisplayName', '参考值');
title('主轴扭矩 - 等效疲劳载荷');
xlabel('风机编号');
ylabel('等效疲劳载荷');
legend('show');
grid on;

% 主轴扭矩的累积疲劳损伤值与参考值对比
subplot(2,1,2);
plot(1:num_wind_turbines, calculated_fatigue_damage_torque, 'b', 'DisplayName', '计算值');
hold on;
plot(1:num_wind_turbines, ref_fatigue_damage_torque, 'r--', 'DisplayName', '参考值');
title('主轴扭矩 - 累积疲劳损伤值');
xlabel('风机编号');
ylabel('累积疲劳损伤值');
legend('show');
grid on;

% ================== 塔架推力的等效疲劳载荷与参考值对比 ==================
figure;
subplot(2,1,1);
plot(1:num_wind_turbines, calculated_equivalent_load_thrust, 'b', 'DisplayName', '计算值');
hold on;
plot(1:num_wind_turbines, ref_equivalent_load_thrust, 'r--', 'DisplayName', '参考值');
title('塔架推力 - 等效疲劳载荷');
xlabel('风机编号');
ylabel('等效疲劳载荷');
legend('show');
grid on;

% 塔架推力的累积疲劳损伤值与参考值对比
subplot(2,1,2);
plot(1:num_wind_turbines, calculated_fatigue_damage_thrust, 'b', 'DisplayName', '计算值');
hold on;
plot(1:num_wind_turbines, ref_fatigue_damage_thrust, 'r--', 'DisplayName', '参考值');
title('塔架推力 - 累积疲劳损伤值');
xlabel('风机编号');
ylabel('累积疲劳损伤值');
legend('show');
grid on;

% ================== 计算平均相对误差 ==================
% 主轴扭矩平均相对误差
mean_relative_error_equivalent_load_torque = mean(abs((calculated_equivalent_load_torque - ref_equivalent_load_torque) ./ ref_equivalent_load_torque));
mean_relative_error_fatigue_damage_torque = mean(abs((calculated_fatigue_damage_torque - ref_fatigue_damage_torque) ./ ref_fatigue_damage_torque));

% 塔架推力平均相对误差
mean_relative_error_equivalent_load_thrust = mean(abs((calculated_equivalent_load_thrust - ref_equivalent_load_thrust) ./ ref_equivalent_load_thrust));
mean_relative_error_fatigue_damage_thrust = mean(abs((calculated_fatigue_damage_thrust - ref_fatigue_damage_thrust) ./ ref_fatigue_damage_thrust));

% 输出平均相对误差
fprintf('主轴扭矩的等效疲劳载荷平均相对误差: %.2f%%\n', mean_relative_error_equivalent_load_torque * 100);
fprintf('主轴扭矩的累积疲劳损伤值平均相对误差: %.2f%%\n', mean_relative_error_fatigue_damage_torque * 100);
fprintf('塔架推力的等效疲劳载荷平均相对误差: %.2f%%\n', mean_relative_error_equivalent_load_thrust * 100);
fprintf('塔架推力的累积疲劳损伤值平均相对误差: %.2f%%\n', mean_relative_error_fatigue_damage_thrust * 100);


function [ranges, means, counts] = three_point_rainflow(load_data, NPeaks)
    % 输入： load_data 是原始的载荷-时间历程数据
    %        NPeaks 是波峰的最大数量
    % 输出： ranges 是识别出的循环幅值
    %        means 是识别出的循环均值
    %        counts 是每个循环的次数，等于循环起点与终点的索引差
    
    % 使用MATLAB内置的findpeaks函数识别波峰和波谷
    [peaks, peak_indices] = findpeaks(load_data, 'NPeaks', NPeaks);       % 限制波峰识别数量
    [valleys, valley_indices] = findpeaks(-load_data, 'NPeaks', NPeaks);  % 限制波谷识别数量（取负）

    valleys = -valleys;  % 取反，恢复波谷的原始值
    
    % 合并波峰和波谷
    extremes = [peaks; valleys];  % 将波峰和波谷的值组合
    indices = [peak_indices; valley_indices];  % 将波峰和波谷的索引组合
    
    % 按索引排序，确保极值按时间顺序排列
    [indices, sort_order] = sort(indices);
    extremes = extremes(sort_order);  % 按相应索引顺序排序极值
    
    % 初始化变量
    ranges = [];
    means = [];
    counts = [];
    
    % 应用三点式雨流计数法识别循环
    while length(extremes) >= 3
        i = 1;
        loop_detected = false;  % 标志是否识别出循环
        
        while i <= length(extremes) - 2
            % 计算 ΔS1 和 ΔS2
            deltaS1 = abs(extremes(i+1) - extremes(i));
            deltaS2 = abs(extremes(i+2) - extremes(i+1));
            
            if deltaS1 <= deltaS2
                % 识别出一个闭合循环，记录循环幅值和均值
                S_ai = deltaS1;  % 循环幅值
                S_mi = (extremes(i) + extremes(i+1)) / 2;  % 循环均值
                
                % 存储幅值、均值和计数（起点和终点的索引差）
                ranges(end+1) = S_ai;
                means(end+1) = S_mi;
                counts(end+1) = abs(indices(i+1) - indices(i))+2;  % 记录索引差
                
                % 移除闭合循环中的两个点（i和i+1），形成新的序列
                extremes(i:i+1) = [];
                indices(i:i+1) = [];
                
                % 标志设置为识别出循环
                loop_detected = true;
                
                % 重置 i 为1，继续从头处理剩余数据
                i = 1;
            else
                % 如果没有闭合循环，移动到下一个点
                i = i + 1;
            end
        end
        
        if ~loop_detected
            total_diff = 0;
            total_mean = 0;
            total_count = 0;
            
            for j = 1:(length(extremes) - 1)
                S_ai = abs(extremes(j+1) - extremes(j));  % 计算幅值
                S_mi = (extremes(j+1) + extremes(j)) / 2;  % 计算均值
                count = abs(indices(j+1) - indices(j));  % 记录索引差

                % 累加幅值、均值和计数
                total_diff = total_diff + S_ai;
                total_mean = total_mean + S_mi;
                total_count = total_count + count;
            end
            
            % 将剩余部分作为一个完整的循环处理，取均值
            ranges(end+1) = total_diff / (length(extremes) - 1);  % 平均幅值
            means(end+1) = total_mean / (length(extremes) - 1);  % 平均均值
            counts(end+1) = total_count;  % 计数为索引差的总和
            
            break;  % 结束循环
        end
    end
end



% =========== Goodman修正公式 ===========
function L_i = goodman_correction(S_ai, S_mi, sigma_b)
    L_i = S_ai / (1 - S_mi / sigma_b);  % 根据Goodman修正公式计算L_i
end

% =========== 计算等效疲劳载荷 ===========
function [L_N, L_i_vector] = calculate_equivalent_load(L_i, counts, m, design_life)
    L_i_vector = (L_i.^m) .* counts;  % 求和前的向量
    L_N = (sum(L_i_vector) / design_life)^(1/m);  % 根据公式计算等效疲劳载荷
end

% =========== 计算累积疲劳损伤 ===========
function [D, D_vector] = calculate_fatigue_damage(L_i,  counts_torque, C, m)
    N_F = round(exp(-m * log(L_i) + log(C)));  % 根据S-N曲线公式计算疲劳寿命 N_F
    D_vector = 1 ./ N_F;  % 求和前的向量
    D = sum(D_vector);  % 累积损伤值
end

三、问题二模型建立与求解

3.1 问题二求解思路
随着全球对可再生能源需求的不断增长，风力发电已经成为满足能源需求的主要可再生能源之一。风力发电机在运行时，受风力、气动载荷以及结构载荷的影响，长时间运转可能导致设备的疲劳损伤和性能下降。为了准确评估风机的运行状态、延长其使用寿命，必须对风机的关键力学参数进行实时监测和估算。塔架推力和主轴扭矩是两个重要的参数，直接关系到风机的结构稳定性和功率输出性能。

通常情况下，风速与风机的发电功率之间存在正相关关系。当风机在理想条件下运行时，风能被风机叶片最大化地捕获并转化为电能。此时，塔架推力和主轴扭矩达到最小值。然而，当风速过高或发电功率不足时，多余的风能作用于风机系统，会导致塔架推力和主轴扭矩的增加，进而加剧风机的疲劳损伤。因此，准确估算风机在不同风速条件下的塔架推力和主轴扭矩，对于风机的设计、运行管理和寿命评估具有重要意义。

本章旨在建立一个基于风速和功率参考值的数学模型，来估算风机的塔架推力和主轴扭矩。我们结合空气动力学理论和能量守恒原理，推导了相应的力学公式，并通过100台风机的监测数据进行验证。本文将详细介绍模型的构建、求解过程以及与实际参考值的误差分析。
3.2 受力分析
风力发电机受力分析涉及风力对叶片的推力、气动扭矩和主轴的扭矩传递。通过调整叶片迎角、叶片半径以及控制转速，可以优化风力发电机的输出功率并保持系统的稳定性。

（1）主轴扭矩分析

在风力发电机运行过程中，叶片旋转带动主轴产生扭矩。这一扭矩由风速和叶片的空气动力学特性决定。图片1中展示了主轴旋转的方向，主轴上的扭矩是由叶片旋转带来的空气动力学作用力产生的，直接作用于主轴并被传递给发电机。

（2）塔架推力分析

塔架推力是风力对叶片产生的作用力，这一作用力通过叶片传递到整个塔架。塔架推力在图片1中展示为叶片前方的推力方向，它是沿着风向作用的。塔架必须承受风产生的推力并保证整个结构的稳定性。

（3）迎角与气动力矩分析

图片1中还标注了叶片的迎角和气动扭矩。迎角是指叶片相对于来流风的角度，它直接影响气动性能和产生的力矩。通过控制叶片的迎角，可以优化气动效率，从而调整主轴扭矩和发电功率。

叶片旋转时，不仅在旋转平面内产生推力，还会因空气动力效应产生气动扭矩。气动扭矩决定了叶片如何受到风力影响，并最终通过传动系统转化为电能。

（4）叶片几何参数的影响

叶片半径是另一个重要参数，它决定了风轮扫掠面积的大小。图片1中标明了叶片半径，较大的叶片半径意味着更大的扫掠面积和更高的捕风能力，从而产生更大的推力和扭矩。

3.3 模型求解与分析

（1）风能与发电功率的关系

风力发电机的工作原理是将风能转化为机械能，然后通过发电机将机械能转化为电能。因此，风机的功率输出直接与风速及其所捕获的风能相关。风的动能公式为：

其中ρ是空气密度,通常取1.225lg/m^3；A是风轮的扫掠面积，等于A=πr^2，其中r为风机叶片的半径；V是风速。风轮从风中捕获的能量与其功率系数C_p有关，风机提取的有效功率可表示为：

其中，C_p是风机的功率系数，取决于风机的设计和运行条件。功率系数的理论最大值为 0.593 (Betz定律)，但实际运行中功率系数通常在 0.35 到 0.5 之间。发电机的输出功率P_ref与风轮提取的功率P_w通过发电机效率η相关联。发电功率参考值P_ref的表达式为：

在本模型中，已知风速V和功率参考值P_ref,我们通过风轮的功率系数C_p以及发电机效率η来反推出塔架推力和主轴扭矩。

（2）主轴扭矩的估算

主轴扭矩是风机从风中捕获的机械能通过叶片传递至主轴的扭矩。主轴扭矩可以通过风机的发电功率与转子转速之间的关系推导出来。风轮上的机械功率P_mech与主轴扭矩T_s和角速度ω_r的关系为：

其中：

⋅Pmech是主轴上的机械功率；

⋅T_s是主轴扭矩；

⋅ ω_r是风轮的角速度。

因此，主轴扭矩可以通过发电功率P_ref和角速度ω_r计算得到：

角速度ω_r可通过风速V和风机叶片的半径r以及尖速比λ计算得到，具体关系为：

其中，λ为风轮的尖速比，是风机叶尖速度与风速之比，通常取值在6到8之间。因此，主轴扭矩的最终计算公式为：

该公式表明，主轴扭矩与功率参考值、叶片半径和风速成正比。

（3）塔架推力的估算

塔架推力是指风机叶片在旋转过程中，风推动叶片产生的水平推力。塔架推力的大小决定了风机塔架所承受的载荷，对风机结构的稳定性具有重要影响。塔架推力可以通过气动力学推导得出，其公式为：

其中·C_t是推力系数，通常取决于风机的设计和运行状况，取值范围一般为 0.7 到 0.9;⋅ρ为空气密度；A=πr^2是风轮的扫掠面积；·V为风速。

推力系数C_t与功率系数C_p之间存在一定的关联，一般来说，当风机功率系数较大时，推力系数也较大。因此，通过风速和功率参考值可以估算出塔架推力的大小。

（4）误差分析

为了评估模型的准确性，本文将模型估算得到的主轴扭矩和塔架推力与风机实际参考值进行对比。通过计算估算值与参考值之间的误差平方和，来衡量模型的精度。误差平方和的计算公式为：

其中：

Tsi和Fti分别为第i台风机的估算主轴扭矩和塔架推力；

·Trefi和Frefi分别为实际参考主轴扭矩和塔架推力。通过分析误差平方和，我们可以确定模型在不同风速和功率工况下的表现，并找出可能导致误差的原因。

3.4 模型求解分析

在该问题中，优化的目标是找到能够最小化主轴扭矩误差和塔架推力误差的最优参数组合。对于每一个风机，利用双目标优化算法计算最优参数，

1. 被优化的参数

功率系数C_p:表示风能转化为机械能的效率，其范围通常为[0.3,0.5]。

空气密度ρ: 空气密度会因风场位置和环境条件而变化，影响推力和扭矩，优化范围为[1.1,1.3]kg/m。

叶片半径R: 叶片的长短直接影响风轮的扫风面积，进而影响推力和扭矩的大小，优化范围为[60,65]米。

2. 优化方法

为了更好地优化风机的功率系数、叶片半径和空气密度三个可调节参数，我们采用了三种多目标优化方法来寻找这些参数的最优解。这三种方法分别是 NSGA-II（非支配排序遗传算法）、加权和法以及 Pareto 前沿构造法。以下对每种优化方法的特点和应用进行详细描述。

（1）NSGA-II（非支配排序遗传算法）

NSGA-II 是一种经典的基于遗传算法的多目标优化方法，能够在不损失多个目标函数的情况下，生成一组 Pareto 最优解集。在风机系统的优化中，我们同时最小化两个目标函数：主轴扭矩误差和塔架推力误差。

NSGA-II 的工作流程：

初始化种群：随机生成若干个初始解，种群中的每个个体代表一组参数组合，包括功率系数、叶片半径和空气密度。

目标函数评估：对于种群中的每个个体，计算两个目标函数值，即主轴扭矩误差和塔架推力误差。误差通过估算值与参考值的平方和计算。

选择、交叉和变异操作：通过基于适应度的选择策略，选择出优质个体进行交叉和变异操作，生成下一代解。

非支配排序和拥挤距离计算：对新生成的个体进行非支配排序，找出表现优异的解；同时通过拥挤距离对解的分布进行调整，确保解在 Pareto 前沿上的均匀分布。

迭代：重复选择、交叉和变异操作，逐渐逼近 Pareto 最优解集，直到达到设定的终止条件（如计算时间达到上限）。

NSGA-II 的优势在于它能够同时优化多个目标，并提供一个 Pareto 最优解集，供决策者根据需要在解集中选择最优解。

（2）加权和法 (Weighted Sum Method)

加权和法是一种将多个目标函数合并为一个目标函数的多目标优化方法。该方法通过引入权重参数 w_1和w_2,将主轴扭矩误差和塔架推力误差的影响合并成一个单一目标函数。优化问题转化为：

其中：

· ETs表示主轴扭矩误差的平方和。·EFt表示塔架推力误差的平方和。· w1和w2是相应的权重，用于调整两个目标之间的相对重要性。

加权和法的工作流程：

设置权重：首先设定权重w1和w2的初始值。权重的调整可以反映不同优化目标的优先级。例如，若更重视最小化主轴扭矩误差，则增大w1的值。

单目标优化：将合并后的目标函数作为单一目标，使用如 fminsearch 这样的单目标优化算法进行求解。通过迭代，找到能够最小化加权和误差的参数组合。

调整权重并重复：若优化结果不符合预期，则可调整权重 w1和w2，重复上述过程，直到得到最优解。

加权和法的优点是实现了从多目标到单一目标的转化，简化了优化过程。然而，该方法对权重的选择较为敏感，若权重设置不当，可能无法找到所有的 Pareto 最优解。

（3）Pareto 前沿构造法（Pareto-front Method）

Pareto 前沿构造法是一种直接寻找 Pareto 最优解集的方法。在多目标优化问题中，Pareto 最优解集代表了所有在不同目标之间无法互相超越的解。这意味着任何一个 Pareto 最优解，在某一目标上变得更优的同时，会导致另一个目标变得更差。

Pareto 前沿构造法的工作流程：

初始化参数范围：设定待优化参数的范围，包括 CpC_pCp、叶片半径 RRR 和空气密度 ρ\rhoρ。

目标函数评估：对于每一组参数组合，计算两个目标函数：主轴扭矩误差和塔架推力误差。

Pareto 排序：对所有候选解进行 Pareto 排序，找出那些在两个目标上均无法被其他解支配的解集。

构造 Pareto 前沿：从所有候选解中构造 Pareto 前沿，生成一组 Pareto 最优解。这些解在不同目标函数之间表现出合理的权衡关系，供决策者选择。

优势：Pareto 前沿构造法能够直接生成一组 Pareto 最优解，允许在不同目标之间进行灵活的权衡选择。这在多目标优化问题中非常有用。

3. 优化结果

下图给出了针对100个风机，三种方法计算得到全部时刻估算值与参照值之差的平方和统计图。表格给出了针对每种算法，100个风机平方和的总和表格，根据表格中的全部时刻估算值与参照值之差的平方和，加权和法表现较好，因此将加权和法获得的全部时刻估算值填入附件6中。

	加权和法	NSGA-II	Pareto-front
所有风机主轴扭矩总误差平方和	2.3317e+22	2.3317e+22	2.3317e+22
所有风机塔架推力总误差平方和	1.6234e+15	1.0655e+16	1.8554e+16

3.5 结果分析

通过对100台风机的数据分析，可以看出，本文建立的数学模型能够有效估算风机的主轴扭矩和塔架推力。模型的误差平方和较小，特别是在额定风速附近，估算结果与实际参考值非常接近。这表明，基于风速和功率参考值的模型在大多数工况下是可靠的。

然而，在风速超过额定值时，风机的控制策略通常会限制叶片的转速以保护设备，这可能导致主轴扭矩的变化不再遵循理论计算结果。因此，未来的模型改进可以考虑引入更多风机的控制策略参数，以提高模型的适应性。

此外，风机运行时受到的载荷不仅与风速和功率有关，还与风切变、湍流强度等气象因素相关。这些因素可能会影响塔架推力和主轴扭矩的实际值。在未来的研究中，可以通过结合更多的气象数据来改进模型，从而进一步提高估算精度。

3.6 参考代码

clc,clear,close all;
% ===================== MATLAB多转子风机控制代码 ====================
% 调用三个优化函数
weighted_sum_optimization();
nsga2_optimization();
pareto_front_optimization();

% ===================== NSGA-II 优化 ====================
function nsga2_optimization()
    % 加载数据
    Pref = readmatrix('sheet1.csv');
    V = readmatrix('sheet2.csv');
    Tshaft_ref = readmatrix('sheet3.csv');
    Ft_ref = readmatrix('sheet4.csv');
    omega_r = readmatrix('sheet6.csv');
    omega_f = readmatrix('sheet7.csv');

    % 固定参数
    generator_efficiency = 0.9; % 发电机效率
    gear_ratio = 100;           % 齿轮箱速比

    % 参数范围定义，只优化Cp、blade_radius和air_density
    lb = [0.3, 60, 1.1]; % 下界：Cp, blade_radius, air_density
    ub = [0.5, 65, 1.3]; % 上界

    % 设定时间限制为5分钟
    options = optimoptions('gamultiobj', 'MaxTime', 300, 'PopulationSize', 100, 'Display', 'iter');

    torque_errors = zeros(100, 1); % 用于存储每个风机的主轴扭矩总误差平方和
    thrust_errors = zeros(100, 1); % 用于存储每个风机的塔架推力总误差平方和

    for turbine = 1:100
        % 目标函数，使用固定的generator_efficiency和gear_ratio
        objective_func_nsga2 = @(x) calculate_torque_and_thrust_fixed(x, turbine, Pref, V, Tshaft_ref, Ft_ref, omega_r, omega_f, generator_efficiency, gear_ratio);

        % 使用 NSGA-II 优化
        [optimal_params, ~] = gamultiobj(objective_func_nsga2, 3, [], [], [], [], lb, ub, options);

        % 计算并保存结果
        [torque_errors(turbine), thrust_errors(turbine)] = save_optimization_results(optimal_params, turbine, Pref, V, Tshaft_ref, Ft_ref, omega_r, omega_f, generator_efficiency, gear_ratio);
    end

    % 绘制100个风机主轴扭矩和塔架推力的总误差平方和
    figure;
    subplot(2, 1, 1);
    bar(torque_errors);
    title('NSGA-II: 主轴扭矩总误差平方和');
    xlabel('风机编号');
    ylabel('总误差平方和');

    subplot(2, 1, 2);
    bar(thrust_errors);
    title('NSGA-II: 塔架推力总误差平方和');
    xlabel('风机编号');
    ylabel('总误差平方和');

    % 计算所有风机的误差平方和总和
    total_torque_error = sum(torque_errors);
    total_thrust_error = sum(thrust_errors);
    fprintf('NSGA-II: 所有风机主轴扭矩总误差平方和: %.4e\n', total_torque_error);
    fprintf('NSGA-II: 所有风机塔架推力总误差平方和: %.4e\n', total_thrust_error);
end

% ===================== 加权和法优化 ====================
function weighted_sum_optimization()
    % 加载数据
    Pref = readmatrix('sheet1.csv');
    V = readmatrix('sheet2.csv');
    Tshaft_ref = readmatrix('sheet3.csv');
    Ft_ref = readmatrix('sheet4.csv');
    omega_r = readmatrix('sheet6.csv');
    omega_f = readmatrix('sheet7.csv');

    % 固定参数
    generator_efficiency = 0.9; % 发电机效率
    gear_ratio = 100;           % 齿轮箱速比

    % 权重
    w1 = 0.5;
    w2 = 0.5;

    torque_errors = zeros(100, 1); % 用于存储每个风机的主轴扭矩总误差平方和
    thrust_errors = zeros(100, 1); % 用于存储每个风机的塔架推力总误差平方和

    for turbine = 1:100
        % 目标函数，只优化Cp, blade_radius, air_density
        objective_func_weighted_sum = @(x) calculate_weighted_error_fixed(x, turbine, Pref, V, Tshaft_ref, Ft_ref, omega_r, omega_f, generator_efficiency, gear_ratio, w1, w2);

        % 初始参数
        init_params = [0.4, 63, 1.225]; % 初始值: Cp, blade_radius, air_density
        
        % 设置优化选项
        options = optimset('MaxIter', 300, 'Display', 'iter');

        % 使用 fminsearch 优化
        optimal_params = fminsearch(objective_func_weighted_sum, init_params, options);

        % 计算并保存结果
        [torque_errors(turbine), thrust_errors(turbine)] = save_optimization_results(optimal_params, turbine, Pref, V, Tshaft_ref, Ft_ref, omega_r, omega_f, generator_efficiency, gear_ratio);
    end

    % 绘制100个风机主轴扭矩和塔架推力的总误差平方和
    figure;
    subplot(2, 1, 1);
    bar(torque_errors);
    title('加权和法: 主轴扭矩总误差平方和');
    xlabel('风机编号');
    ylabel('总误差平方和');

    subplot(2, 1, 2);
    bar(thrust_errors);
    title('加权和法: 塔架推力总误差平方和');
    xlabel('风机编号');
    ylabel('总误差平方和');

    % 计算所有风机的误差平方和总和
    total_torque_error = sum(torque_errors);
    total_thrust_error = sum(thrust_errors);
    fprintf('加权和法: 所有风机主轴扭矩总误差平方和: %.4e\n', total_torque_error);
    fprintf('加权和法: 所有风机塔架推力总误差平方和: %.4e\n', total_thrust_error);
end

% ===================== Pareto-front 优化 ====================
function pareto_front_optimization()
    % 加载数据
    Pref = readmatrix('sheet1.csv');
    V = readmatrix('sheet2.csv');
    Tshaft_ref = readmatrix('sheet3.csv');
    Ft_ref = readmatrix('sheet4.csv');
    omega_r = readmatrix('sheet6.csv');
    omega_f = readmatrix('sheet7.csv');

    % 固定参数
    generator_efficiency = 0.9; % 发电机效率
    gear_ratio = 100;           % 齿轮箱速比

    % 参数范围定义，只优化Cp, blade_radius和air_density
    lb = [0.3, 60, 1.1]; % 下界
    ub = [0.5, 65, 1.3]; % 上界

    % 设置时间限制为5分钟
    options = optimoptions('paretosearch', 'MaxTime', 300, 'Display', 'iter');

    torque_errors = zeros(100, 1); % 用于存储每个风机的主轴扭矩总误差平方和
    thrust_errors = zeros(100, 1); % 用于存储每个风机的塔架推力总误差平方和

    for turbine = 1:100
        % 目标函数
        objective_func_pareto = @(x) calculate_torque_and_thrust_fixed(x, turbine, Pref, V, Tshaft_ref, Ft_ref, omega_r, omega_f, generator_efficiency, gear_ratio);

        % 使用 Pareto-front 优化
        [optimal_params, ~] = paretosearch(objective_func_pareto, 3, [], [], [], [], lb, ub, [], options);

        % 计算并保存结果
        [torque_errors(turbine), thrust_errors(turbine)] = save_optimization_results(optimal_params, turbine, Pref, V, Tshaft_ref, Ft_ref, omega_r, omega_f, generator_efficiency, gear_ratio);
    end

    % 绘制100个风机主轴扭矩和塔架推力的总误差平方和
    figure;
    subplot(2, 1, 1);
    bar(torque_errors);
    title('Pareto-front: 主轴扭矩总误差平方和');
    xlabel('风机编号');
    ylabel('总误差平方和');

    subplot(2, 1, 2);
    bar(thrust_errors);
    title('Pareto-front: 塔架推力总误差平方和');
    xlabel('风机编号');
    ylabel('总误差平方和');

    % 计算所有风机的误差平方和总和
    total_torque_error = sum(torque_errors);
    total_thrust_error = sum(thrust_errors);
    fprintf('Pareto-front: 所有风机主轴扭矩总误差平方和: %.4e\n', total_torque_error);
    fprintf('Pareto-front: 所有风机塔架推力总误差平方和: %.4e\n', total_thrust_error);
end

% ===================== 计算与保存结果的函数 ====================
function [torque_error_sum_of_squares, thrust_error_sum_of_squares] = save_optimization_results(optimal_params, turbine, Pref, V, Tshaft_ref, Ft_ref, omega_r, omega_f, generator_efficiency, gear_ratio)
    % 计算结果并与参考值对比
    num_time_steps = size(Pref, 1);
    Cp = optimal_params(1);
    blade_radius = optimal_params(2);
    air_density = optimal_params(3);

    Tshaft_est = zeros(num_time_steps, 1);
    Ft_est = zeros(num_time_steps, 1);

    for t = 1:num_time_steps
        P_ref = Pref(t, turbine);
        wind_speed = V(t, turbine);
        omega_r_turbine = omega_r(t, turbine);
        omega_f_turbine = omega_f(t, turbine);

        % 主轴扭矩计算
        Tshaft_est(t) = (P_ref / (generator_efficiency * omega_r_turbine)) * gear_ratio;
        Tshaft_est(t) = min(max(Tshaft_est(t), 0), 10^9);

        % 塔架推力计算
        swept_area = pi * blade_radius^2;
        Ft_est(t) = 0.5 * air_density * swept_area * Cp * (wind_speed^2);
        Ft_est(t) = min(max(Ft_est(t), 0), 10^7);
    end

    % 计算总误差平方和
    torque_error_sum_of_squares = sum((Tshaft_est - Tshaft_ref(:, turbine)).^2);
    thrust_error_sum_of_squares = sum((Ft_est - Ft_ref(:, turbine)).^2);

    % 保存数据到附件6的表格中
    % result_data = [Tshaft_est, Ft_est];
    % writematrix(result_data, '附件6.xlsx', 'Sheet', sprintf('WindTurbine_%d', turbine));
end

% ===================== 主轴扭矩和塔架推力估算的函数 ====================
function [torque_error, thrust_error] = calculate_torque_and_thrust_fixed(params, turbine, Pref, V, Tshaft_ref, Ft_ref, omega_r, omega_f, generator_efficiency, gear_ratio)
    % 提取参数
    Cp = params(1);                   % 功率系数
    blade_radius = params(2);          % 叶片半径
    air_density = params(3);           % 空气密度

    % 初始化估算值
    num_time_steps = size(Pref, 1);
    Tshaft_est = zeros(num_time_steps, 1);
    Ft_est = zeros(num_time_steps, 1);

    for t = 1:num_time_steps
        P_ref = Pref(t, turbine);
        wind_speed = V(t, turbine);
        omega_r_turbine = omega_r(t, turbine);
        omega_f_turbine = omega_f(t, turbine);

        % 主轴扭矩计算
        Tshaft_est(t) = (P_ref / (generator_efficiency * omega_r_turbine)) * gear_ratio;
        Tshaft_est(t) = min(max(Tshaft_est(t), 0), 10^9);  % 限制估算值

        % 塔架推力计算
        swept_area = pi * blade_radius^2;
        Ft_est(t) = 0.5 * air_density * swept_area * Cp * (wind_speed^2);
        Ft_est(t) = min(max(Ft_est(t), 0), 10^7);  % 限制推力估算值
    end

    % 计算主轴扭矩和塔架推力的误差平方和
    torque_error = sum((Tshaft_est - Tshaft_ref(:, turbine)).^2);
    thrust_error = sum((Ft_est - Ft_ref(:, turbine)).^2);
end

function weighted_error = calculate_weighted_error_fixed(params, turbine, Pref, V, Tshaft_ref, Ft_ref, omega_r, omega_f, w1, w2, generator_efficiency, gear_ratio)
    % 提取可调参数
    Cp = params(1);                   % 功率系数
    blade_radius = params(2);          % 风机叶片半径
    air_density = params(3);           % 空气密度

    % 初始化估算值
    num_time_steps = size(Pref, 1);
    Tshaft_est = zeros(num_time_steps, 1);  % 估算的主轴扭矩
    Ft_est = zeros(num_time_steps, 1);      % 估算的塔架推力

    % 循环计算每个时间步的主轴扭矩和塔架推力
    for t = 1:num_time_steps
        P_ref = Pref(t, turbine);          % 功率调度指令
        wind_speed = V(t, turbine);        % 输入风速
        omega_r_turbine = omega_r(t, turbine);  % 低速轴转速
        omega_f_turbine = omega_f(t, turbine);  % 高速轴转速

        % ======== 主轴扭矩估算 ========
        if wind_speed > 0
            % 基于功率、发电机效率和齿轮箱速比估算主轴扭矩
            Tshaft_est(t) = (P_ref / (generator_efficiency * omega_r_turbine)) * gear_ratio;
            Tshaft_est(t) = min(max(Tshaft_est(t), 0), 10^9);  % 限制扭矩估算值在合理范围内
        else
            Tshaft_est(t) = 0;
        end

        % ======== 塔架推力估算 ========
        swept_area = pi * blade_radius^2;  % 计算叶片扫过的面积
        Ft_est(t) = 0.5 * air_density * swept_area * Cp * (wind_speed^2);  % 基于风速计算推力
        Ft_est(t) = min(max(Ft_est(t), 0), 10^7);  % 限制推力估算值在合理范围内
    end

    % 计算主轴扭矩和塔架推力的误差平方和
    torque_error = sum((Tshaft_est - Tshaft_ref(:, turbine)).^2);
    thrust_error = sum((Ft_est - Ft_ref(:, turbine)).^2);

    % 加权总误差（加权和法）
    weighted_error = w1 * torque_error + w2 * thrust_error;
end


clc,clear,close all;
% ===================== MATLAB多转子风机控制代码 ====================
% 调用三个优化函数
weighted_sum_optimization();




% ===================== 加权和法优化 ====================
function weighted_sum_optimization()
    % 加载数据
    Pref = readmatrix('sheet1.csv');
    V = readmatrix('sheet2.csv');
    Tshaft_ref = readmatrix('sheet3.csv');
    Ft_ref = readmatrix('sheet4.csv');
    omega_r = readmatrix('sheet6.csv');
    omega_f = readmatrix('sheet7.csv');

    % 权重
    w1 = 0.5;
    w2 = 0.5;

    torque_errors = zeros(100, 1); % 用于存储每个风机的主轴扭矩总误差平方和
    thrust_errors = zeros(100, 1); % 用于存储每个风机的塔架推力总误差平方和
    result1=zeros(2000,100);
    result2=zeros(2000,100);
    for turbine = 1:100
        % 目标函数
        objective_func_weighted_sum = @(x) calculate_weighted_error(x, turbine, Pref, V, Tshaft_ref, Ft_ref, omega_r, omega_f, w1, w2);

        % 初始参数
        init_params = [0.4, 63, 1.225, 0.9, 100];
        
        % 设置优化选项
        options = optimset('MaxIter', 300, 'Display', 'iter');

        % 使用 fminsearch 优化
        optimal_params = fminsearch(objective_func_weighted_sum, init_params, options);

        % 计算并保存结果
        [result1(:,turbine), result2(:,turbine), torque_errors(turbine), thrust_errors(turbine)] = save_optimization_results(optimal_params, turbine, Pref, V, Tshaft_ref, Ft_ref, omega_r, omega_f);
    
    end
    % 保存数据到附件6的表格中
    %result_data = [Tshaft_est, Ft_est];
    writematrix(result1, '附件6.xlsx', 'Sheet', sprintf('WindTurbine_%d', 1));
    writematrix(result2, '附件6.xlsx', 'Sheet', sprintf('WindTurbine_%d', 2));



    % 绘制100个风机主轴扭矩和塔架推力的总误差平方和
    figure;
    subplot(2, 1, 1);
    bar(torque_errors);
    title('加权和法: 主轴扭矩总误差平方和');
    xlabel('风机编号');
    ylabel('总误差平方和');

    subplot(2, 1, 2);
    bar(thrust_errors);
    title('加权和法: 塔架推力总误差平方和');
    xlabel('风机编号');
    ylabel('总误差平方和');

    % 计算所有风机的误差平方和总和
    total_torque_error = sum(torque_errors);
    total_thrust_error = sum(thrust_errors);
    fprintf('加权和法: 所有风机主轴扭矩总误差平方和: %.4e\n', total_torque_error);
    fprintf('加权和法: 所有风机塔架推力总误差平方和: %.4e\n', total_thrust_error);
end


% ===================== 计算与保存结果的函数 ====================
function [Tshaft_est, Ft_est,torque_error_sum_of_squares, thrust_error_sum_of_squares] = save_optimization_results(optimal_params, turbine, Pref, V, Tshaft_ref, Ft_ref, omega_r, omega_f)
    % 计算结果并与参考值对
    % 计算结果并与参考值对比
    num_time_steps = size(Pref, 1);
    Cp = optimal_params(1);
    blade_radius = optimal_params(2);
    air_density = optimal_params(3);
    generator_efficiency = optimal_params(4);
    gear_ratio = optimal_params(5);

    Tshaft_est = zeros(num_time_steps, 1);
    Ft_est = zeros(num_time_steps, 1);

    for t = 1:num_time_steps
        P_ref = Pref(t, turbine);
        wind_speed = V(t, turbine);
        omega_r_turbine = omega_r(t, turbine);
        omega_f_turbine = omega_f(t, turbine);

        % 主轴扭矩计算
        Tshaft_est(t) = (P_ref / (generator_efficiency * omega_r_turbine)) * gear_ratio;
        Tshaft_est(t) = min(max(Tshaft_est(t), 0), 10^7);

        % 塔架推力计算
        swept_area = pi * blade_radius^2;
        Ft_est(t) = 0.5 * air_density * swept_area * Cp * (wind_speed^2);
        Ft_est(t) = min(max(Ft_est(t), 0), 10^7);
    end

    % 计算总误差平方和
    torque_error_sum_of_squares = sum((Tshaft_est - Tshaft_ref(:, turbine)).^2);
    thrust_error_sum_of_squares = sum((Ft_est - Ft_ref(:, turbine)).^2);

    % 打印结果
    fprintf('风机 %d 主轴扭矩总误差平方和: %.4e\n', turbine, torque_error_sum_of_squares);
    fprintf('风机 %d 塔架推力总误差平方和: %.4e\n', turbine, thrust_error_sum_of_squares);

   end

function weighted_error = calculate_weighted_error(params, turbine, Pref, V, Tshaft_ref, Ft_ref, omega_r, omega_f, w1, w2)
    % 提取输入参数
    Cp = params(1);                  % 功率系数
    blade_radius = params(2);         % 风机叶片半径
    air_density = params(3);          % 空气密度
    generator_efficiency = params(4); % 发电机效率
    gear_ratio = params(5);           % 齿轮箱速比

    % 初始化估算值
    num_time_steps = size(Pref, 1);
    Tshaft_est = zeros(num_time_steps, 1);  % 估算的主轴扭矩
    Ft_est = zeros(num_time_steps, 1);      % 估算的塔架推力

    % 循环计算每个时间步的主轴扭矩和塔架推力
    for t = 1:num_time_steps
        P_ref = Pref(t, turbine);          % 功率调度指令
        wind_speed = V(t, turbine);        % 输入风速
        omega_r_turbine = omega_r(t, turbine);  % 低速轴转速
        omega_f_turbine = omega_f(t, turbine);  % 高速轴转速

        % ======== 主轴扭矩估算 ========
        if wind_speed > 0
            % 基于功率、发电机效率和齿轮箱速比估算主轴扭矩
            Tshaft_est(t) = (P_ref / (generator_efficiency * omega_r_turbine)) * gear_ratio;
            Tshaft_est(t) = min(max(Tshaft_est(t), 0), 10^9);  % 限制扭矩估算值在合理范围内
        else
            Tshaft_est(t) = 0;
        end

        % ======== 塔架推力估算 ========
        swept_area = pi * blade_radius^2;  % 计算叶片扫过的面积
        Ft_est(t) = 0.5 * air_density * swept_area * Cp * (wind_speed^2);  % 基于风速计算推力
        Ft_est(t) = min(max(Ft_est(t), 0), 10^7);  % 限制推力估算值在合理范围内
    end

    % 计算主轴扭矩和塔架推力的误差
    torque_error = sum((Tshaft_est - Tshaft_ref(:, turbine)).^2);
    thrust_error = sum((Ft_est - Ft_ref(:, turbine)).^2);

    % 加权总误差（加权和法）
    weighted_error = w1 * torque_error + w2 * thrust_error;
end

...

参考文献

汤国生.基于SWOT分析法的大学生数学建模创新实践基地建设探索——以江苏科技大学为例[J].高教学刊,2023,9(11):53-56.DOI:10.19980/j.CN23-1593/G4.2023.11.013.
谢雨婧,韩惠丽.北师大版初中数学教材中数学建模的多维度分析[J].教学与管理,2023(09):73-76.
刘志梅.数学建模与高职数学教学的深度融合[J].佳木斯职业学院学报,2023,39(03):152-154.
许亚桃,吴立宝.基于Delphi-AHP高中数学建模教学评价指标体系的研究[J].内江师范学院学报,2023,38(02):113-119.DOI:10.13603/j.cnki.51-1621/z.2023.02.018.
杨本朝,石雅男,段乾恒,李光松,于刚.大学生数学建模竞赛开展全周期教学实践探究[J].大学教育,2023(04):44-46.
黄健,徐斌艳.国际视野下数学建模教与学研究的发展趋势——基于第14届国际数学教育大会的分析[J].数学教育学报,2023,32(01):93-98.

...