基于理论分析与数值模拟结果提出更简洁的旋筒帆升阻力系数回归公式

根据文献阅读与分析，分析旋筒帆的升力系数以及阻力系数随着自旋比、长径比以及端板直径比提出适合进行多元非线性回归的更简洁的升阻力系数表达式，对应的自变量依旧是旋筒帆的长径比、自旋比以及端板直径比，由于旋筒帆的端板直径比过大会导致升阻力系数的非规律性变化，且对应的自旋功率消耗也会急剧增大，违背了节能减排的初衷，因此将端板直径比的范围限制在1.5-3.0之间，并且对旋筒帆的长径比也限制在5.5-7.5之间，自旋比限制在2.0-4.5之间，并对其进行多元非线性回归。提出的应用于多元非线性回归的升阻力系数计算公式如下：

升阻力系数回归公式表达式

队员非线性回归的样本数据如下：

多元非线性回归的样本数据集

长径比	自旋比	端板直径比	升力系数	阻力系数
8.000	0.938	1.500	1.552	0.462
8.000	2.000	1.500	5.379	1.038
8.000	3.012	1.500	8.483	1.865
8.000	4.000	1.500	9.621	2.519
8.000	5.012	1.500	9.517	2.558
8.000	0.938	2.000	1.862	0.404
8.000	2.000	2.000	6.000	0.962
8.000	3.012	2.000	9.000	1.769
8.000	4.000	2.000	10.241	2.385
8.000	5.012	2.000	10.552	2.385
8.000	0.938	3.000	2.069	0.385
8.000	2.000	3.000	6.517	0.827
8.000	3.012	3.000	9.828	1.538
8.000	4.000	3.000	11.379	2.058
8.000	5.012	3.000	11.897	2.000
4.000	1.000	2.000	1.733	0.662
4.000	2.000	2.000	5.600	1.457
4.000	3.000	2.000	8.267	2.596
4.000	4.000	2.000	9.467	3.444
4.000	5.000	2.000	9.400	3.815
6.000	1.000	2.000	1.800	0.397
6.000	2.000	2.000	5.600	1.139
6.000	3.000	2.000	8.400	2.040
6.000	4.000	2.000	9.600	2.808
6.000	5.000	2.000	9.667	2.967
8.000	1.000	2.000	1.800	0.397
8.000	2.000	2.000	5.600	0.954
8.000	3.000	2.000	8.533	1.748
8.000	4.000	2.000	9.867	2.358
8.000	5.000	2.000	9.933	2.384

采用控制回归系数上限限的方法求得到最佳的上下限系数，再基于matlab自带的lsqcurvefit函数，进行多元非线性回归，可以得到如下计算结果：

（1）升力系数回归计算结果：

（2）阻力系数回归计算结果:

最终非线性回归得到的公式各项系数结果为：

系数	升力系数回归模型系数值	阻力系数回归模型系数值
b(1)	0.0457	-0.3562
b(2)	-0.5673	2.7702
b(3)	1.9999	-2.6402
b(4)	-0.1736	2.6374
b(5)	0.02629	-0.4795
b(6)	1.62986	7.1014
b(7)	-0.99999	4.9987
b(8)	-0.31382	-6.4054

可见新提出的非线性回归公式的回归效果很好，误差集中在5%-10%左右，满足旋筒帆的升阻力系数的快速预报任务要求精度。（当然不同精度下回归公式结果不同）

下面将对比非样本数据下的数值模拟和实验的结果与多项式回归计算结果（2024，Optimal Design of Rotor Sails Based on Environmental Conditions and Cost，Cem Guzelbulut *等人）以及新公式快速预报进行对比，结果如下：

升力系数预报结果对比图

阻力系数预报结果对比图

可见预报结果较好，根据回归分析得到的各项系数可知，新的多元非线性回归公式计算不会出现负值的情况，能够更好的预测不同自旋比、长径比和端板直径比下的旋筒帆升阻力系数值，并应用与船舶的旋筒帆的参数设计以及后续的航线规划中，提供更准确、可靠的参考依据。

% 创建一个矩阵来存储数据
data = xlsread('222.xlsx', 1, 'A1:D66');

% 将数据分成输入和输出变量
X = data(:, 1:3); % 长径比、自旋比、端板直径比
Y = data(:, 4);   % 升力系数

% 设置交叉验证的折数
numFolds = 5;
cv = cvpartition(size(X,1), 'KFold', numFolds);

% 初始化最佳参数和性能
bestRsq = -0.05;
bestLb = [];
bestUb = [];

% 生成可能的上界和下界的参数范围
lb_range = -1:0.05:1; % 设置下界范围
ub_range = 1:0.05:2; % 设置上界范围

% 通过交叉验证找到最佳的上界和下界
for lb = lb_range
    for ub = ub_range
        fun = @(b, x) (b(1)*x(:,2).^4 + b(2)*x(:,2).^3 + b(3)*x(:,2).^2 + b(4)*x(:,2).^1)...
            .*(1.05).^(b(5).*x(:,1) + b(6).*x(:,1).^0.5).*(0.95).^(b(7).*x(:,3) + b(8).*x(:,3).^2);
        % 使用交叉验证评估模型性能
        rsq_values = zeros(numFolds, 1);
        for fold = 1:numFolds
            trainIdx = cv.training(fold);
            testIdx = cv.test(fold);
            Xtrain = X(trainIdx, :);
            Ytrain = Y(trainIdx);
            Xtest = X(testIdx, :);
            Ytest = Y(testIdx);

            % 使用 lsqcurvefit 进行拟合
            beta0 = zeros(8, 1);
            lb_vec = ones(8, 1) * lb;
            ub_vec = ones(8, 1) * ub;
            
            % 添加惩罚项
            fun_with_penalty = @(b, x)fun(b, x) + penalty(b);
            
            beta = lsqcurvefit(fun_with_penalty, beta0, Xtrain, Ytrain, lb_vec, ub_vec);

            % 计算 R 平方值
            YFit = fun(beta, Xtest);
            residuals = Ytest - YFit;
            SSresid = sum(residuals.^2);
            SStotal = (length(Ytest) - 1) * var(Ytest);
            rsq_values(fold) = 1 - SSresid/SStotal;
        end

        % 取平均 R 平方值
        avgRsq = mean(rsq_values);

        % 更新最佳参数和性能
        if avgRsq > bestRsq
            bestRsq = avgRsq;
            bestLb = lb;
            bestUb = ub;
        end
    end
end

% 输出最佳参数和性能
disp(['最佳下界: ', num2str(bestLb)]);
disp(['最佳上界: ', num2str(bestUb)]);
disp(['最佳平均R平方值: ', num2str(bestRsq)]);

% 使用最佳参数进行拟合
fun =  @(b, x) (b(1)*x(:,2).^4 + b(2)*x(:,2).^3 + b(3)*x(:,2).^2 + b(4)*x(:,2).^1)...
            .*(1.05).^(b(5).*x(:,1) + b(6).*x(:,1).^0.5).*(0.95).^(b(7).*x(:,3) + b(8).*x(:,3).^2);

beta0 = zeros(8, 1);
lb_vec = ones(8, 1) .* bestLb;
ub_vec = ones(8, 1) .* bestUb;

beta = lsqcurvefit(fun, beta0, X, Y, lb_vec, ub_vec);

% 计算拟合值
YFit = fun(beta, X);

% 计算拟合误差
residuals = Y - YFit;
SSresid = sum(residuals.^2);
SStotal = (length(Y) - 1) * var(Y);
rsq = 1 - SSresid/SStotal;
rmse = sqrt(mean(residuals.^2));

% 打印结果
disp('拟合系数:');
disp(beta);
disp(['拟合R平方值: ', num2str(rsq)]);

% 绘制拟合图
figure;
plot(Y, 'ko', 'DisplayName', '原始数据');
hold on;
plot(YFit, 'r-', 'LineWidth', 2, 'DisplayName', '拟合曲线');
xlabel('样本值');
ylabel('拟合值');
title('多元非线性回归拟合结果');
legend('Location', 'best');
grid on;

% 绘制实际值与预测值对比图
figure;
scatter(Y, YFit, 'filled');
hold on;
plot([min(Y), max(Y)], [min(Y), max(Y)], '--k');
xlabel('实际升力系数');
ylabel('预测升力系数');
title('实际升力系数 vs. 预测升力系数');
grid on;
hold off;

% 绘制误差图
figure;
plot(residuals, 'o');
xlabel('样本');
ylabel('拟合误差');
title('拟合误差分析');
grid on;

% 绘制误差分布图
figure;
histogram(residuals,8, 'Normalization', 'probability');
xlabel('预测误差');
ylabel('概率');
title(['预测误差分布 (RMSE = ' num2str(rmse) ')']);
grid on;