2024最新:​鹅优化算法(GO),原理公式详解,附matlab代码

本文介绍了一种新颖的优化算法——鹅优化算法(GOOSE),灵感来源于鹅的休息和觅食行为。该算法通过动态调整搜索策略,具有快速收敛和高精度的特点,已在知名SCI期刊EvolvingSystems上发表。

鹅优化算法(GOOSE AlgorithmGO)是一种群智能优化算法,Rebwar Khalid Hamad提出了一种基于鹅的休息和觅食行为的新型元启发式算法。鹅用一条腿站立,保持平衡,以保护和保护群中的其他个体。该算法通过自适应地调整搜索空间的分辨率和搜索速度,以快速而准确地找到最优解,具有收敛速度快、求解精度高等特点,是一个出色的优化算法。该成果于2024发表在知名SCI期刊JCQ3区Evolving Systems上。

fa4399460ba1140fe7dd0d6866a1f014.png

GO从鹅的休息和觅食行为获得灵感,当鹅听到任何奇怪的声音或动作时,它们会发出响亮的声音来唤醒群中的个体,并保证它们的安全。

鹅群在它们之间创造了一种有吸引力的保护氛围。鹅的动作特点分为三类:1. 在它们休息的时候,大雁成群地聚集在一起,其中一只用单腿保持平衡。2. 偶尔,他抬起一条腿,扛起一块小石头,这样,当他睡着时,石头又掉了下来,鹅就会醒过来。3. 当鹅注意到任何意外的噪音或活动时,它会发出响亮的叫声来警告群中的其他动物,以保证它们的安全。

be7e8d895677e300ebb767e6a647a29f.png

鹅成群地聚集在它们的避难所和休息区。在种群中,会有一只鹅被指派守卫,这只鹅开始通过站立和单腿平衡来执行他的命令。每当鹅睡着时,它的腿或它携带的石头就会掉到地上。这时,岩石的声音传到了群中的另一只鹅。它一听到这个声音就会处于被剥削的状态。因此,其他个体需要一些时间才能听到守护鹅的叫声。

61799b8a95aa6bb62452f59eb6d82ccb.png

1、算法原理

(1)开发阶段

我们首先需要估计鹅脚里的石头的重量,在525公斤之间。通过下式,我们可以在任意迭代中随机求出石头的重量。这个变量指示迭代的次数。

然后,在下式中,我们应该找到当石头落下时到达地球所需的时间time_of_arrive_object。它是随机的,在1到循环中每次迭代的维数之间。

在下式中,我们找到了物体撞击地面并发出声音并传递给群中每只鹅的时间Time_of_Arrive_Soundit

在下一个等式中,找出声音在整个迭代过程中传播并到达群中单个鹅所需的总时间。如该式所示,维度除以总时间量。

为了得到所需的平均时间,我们将总时间除以2。下式解释了步骤。

正如我们之前所讨论的,有一个随机变量rnd负责开发和探索阶段的分布。变量pro的值在[0,1]范围内随机选取。假设变量pro的值大于0.2,且Stone_Weightit大于等于12。在上式中,time_of_arrive_object乘以stone_weight的平方根除以物体的加速度。为了保护和唤醒群体中的个体,这些方程式应该被计算出来。

(2)探索阶段

一旦其中一只鹅醒来,它们就开始尖叫,以保护鸟群中的所有个体,此时开始进入探索阶段。如果变量rnd的值小于0.5,那么应用以下方程。

通过使用randn(1, dim)确保鹅随机探索搜索空间中的其他个体。然而,Minimum_Timealpha变量都被用来提高GOOSE的可搜索性。在下式中,将时间和alpha的最小值乘以一个随机数,然后加到搜索空间中的最佳位置。

Міпітит

其中dim是问题维度的数量,Best_pos是我们在搜索区域中找到的最佳位置。

总而言之,鹅算法随机开始生成种群。然后,在第一次迭代中,它检查种群中的值,以恢复边界外的值。同时,实现目标函数来确定搜索边界内的最佳分数和最佳位置。为了控制开发和勘探阶段,我们使用随机选择值的随机变量rand。如果rand的值大于或等于0.5,则激活探索阶段。

2、结果展示

18e3a3772315b500ba934f9124c2321c.png

844163ba09b720957cfde0a23cb68350.png

70bacaa8055a39df78e88ce4593be288.png

1c1f4b1416dfdd86dfd32a6f6cf34c2d.png

8e9b99655907f1d1d50a4ebed224fa3b.png

可以看到,该算法基本不会陷入局部最优解。搜索效率也还不错!

3、MATLAB核心代码

%% 淘个代码 %%
% 微信公众号搜索:淘个代码,获取更多代码
% 鹅优化算法(GO)
function [Best_score,Best_pos,Cong_Curve]=GOOSE(SearchAgents_no,Max_IT,lb,ub,dim,fobj)
Best_pos=zeros(1,dim);
Best_score=inf;                          %change this to -inf for maximization problems
M_T=inf;
Cong_Curve=zeros(1,Max_IT);
%Initialize the positions of search agents
X=initialization(SearchAgents_no,dim,ub,lb);
Distance_Goose=zeros(SearchAgents_no,dim);
loop=0;                                               % Loop counter
% Main loop
while loop<Max_IT
    for i=1:size(X,1)
        % Return back the search agents that go beyond the boundaries of the search space
        Flag4ub=X(i,:)>ub;
        Flag4lb=X(i,:)<lb;
        X(i,:)=(X(i,:).*(~(Flag4ub+Flag4lb)))+ub.*Flag4ub+lb.*Flag4lb;
        % Calculate objective function for each search agent
        fitness=fobj(X(i,:));
        if fitness<Best_score
            Best_score=fitness;
            Best_pos=X(i,:);
        end
    end
    for i=1:size(X,1)
        pro=rand;
        rnd=rand;
        coe=rand();
        if(coe<=0.17)
            coe;
        else
            coe=0.17;
        end
S_W(i,:)=randi([5,25],1,1);                                         %Eq.(3.1)
        % Time of Arrive object to earth.
T_o_A_O(i,:)=rand(1,dim);                                           %Eq.(3.2)
        % Time of Arrive Sound to Groups
T_o_A_S(i,:)=rand(1,dim);                                           %Eq.(3.3)
T_T=sum(T_o_A_S(i,:))/dim;                                          %Eq.(3.4)
        % Determine Time Average
T_A=T_T/2;                                                          %Eq.(3.5)
        if rnd>=0.5
            if pro>0.2
                if S_W>=12
                    %  Calculate Free Fall Speed
                    F_F_S=T_o_A_O(i,:) *(sqrt(S_W(i,:))/ 9.81);       %Eq.(3.6)
                    S_S=343.2;
                    D_S_T(i,:)=S_S* T_o_A_S(i,:);                     %Eq.(3.7)
                    D_G(i,:)=0.5* D_S_T(i,:);                         %Eq.(3.8)
                    X(i,:)=F_F_S + D_G(i,:)* T_A^2;                     %Eq.(3.9)
                elseif S_W<12
                elseif pro<=0.2
                    %  Calculate Free Fall Speed
                    F_F_S=T_o_A_O(i,:) *(S_W(i,:)/ 9.81);             %Eq.(3.10)
                    S_S=343.2;
                    D_S_T(i,:)=S_S* T_o_A_S(i,:);                     %Eq.(3.7)
                    D_G(i,:)=0.5* D_S_T(i,:);                         %Eq.(3.8)
                    X(i,:)=F_F_S.*D_G(i,:)* T_A^2*coe;                  %Eq.(3.11)
                end
            end
        else
            if M_T>T_T
                M_T=T_T;
            end
            alpha=(2-(loop/(Max_IT/2)));                         %Eq.(3.12)
            %random an awakening
            % exploring wakeup without holding stone.
            X(i,:)=randn(1,dim).*(M_T*alpha)+Best_pos;           %Eq.(3.13)
        end
    end
    loop=loop+1;
    Cong_Curve(loop)=Best_score;
end

参考文献

[1] Hamad R K, Rashid T A. GOOSE algorithm: a powerful optimization tool for real-world engineering challenges and beyond[J]. Evolving Systems, 2024: 1-26.

完整代码获取方式:后台回复关键字:

TGDM199

内容概要:本文档详细介绍了一个基于MATLAB的项目实例,该项目旨在实现一个基于GOOSE-LightGBM的多特征输入单输出数据回归预测系统。通过结合GOOSE(Generalized Optimal One-Versus-One Strategy)和LightGBM算法,项目不仅提升了回归预测的准确性和效率,还能处理大规模多特征数据集。项目包括六个主要阶段:环境准备、算法设计、模型构建、性能评估、GUI设计以及防止过拟合措施,并提供了详细的代码示例。项目特点包括高效的分布式训练、可视化的预测结果分析和全面的数据处理模块,确保模型具备较强的泛化能力与应用范围。 适合人群:具有MATLAB编程经验的研发人员,尤其是从事数据分析与机器学习工作的专业人士,以及希望通过机器学习方法提升预测模型性能的研究人员。 使用场景及目标:①需要处理多维输入特征并在单一目标上进行精确预测的应用;②能源管理、金融市场、医疗诊断、气象预报及交通管理等多个领域的数据驱动型任务,特别是那些需要高频次预测和支持大规模数据集的应用;③旨在优化现有模型性能,提高运算速度及降低成本的企业和机构。 其他说明:在实际应用中,该项目强调了几个关键点,首先是数据的质量管理和预处理,这是确保预测模型准确的基础。其次是模型的稳定性与可靠性,项目通过引入分布式计算、GPU加速、自动调参等技术来提升训练效率并保证长时间内的稳定运行。最后,在安全性和法律合规性方面,项目特别关注用户隐私保护和技术的合法性应用,以满足现代社会严格的监管要求。通过这些努力,项目不仅能解决复杂的回归预测难题,也为后续的技术演进奠定了坚实的基础。
<think> 我们已知输入信号u1, u2和输出信号y,需要确定状态空间模型的最佳阶数。 常用的方法是使用信息准则(如AIC、BIC)或利用子空间方法的奇异值分解(SVD)来确定阶数。 在MATLAB中,我们可以使用系统辨识工具箱(System Identification Toolbox)中的函数,例如`ssest`或`n4sid`,它们可以自动估计状态空间模型并帮助我们选择阶数。 另一种方法是使用`arx`结构来估计模型,然后通过计算不同阶数的损失函数(预测误差的方差)并应用信息准则来选择最佳阶数。 然而,由于我们有多个输入(u1, u2),我们考虑MISO(多输入单输出)系统。我们可以使用状态空间模型,其阶数表示状态向量的维数。 这里提供两种方法: 方法1:使用`ssest`函数,设置阶数范围,然后比较不同阶数下的拟合度(比如AIC或BIC准则)来选择最佳阶数。 方法2:使用子空间方法(如`n4sid`)生成一个模型数组,然后比较这些模型的信息准则值。 我们将使用方法2,因为`n4sid`方法在估计状态空间模型时,可以通过奇异值分解(SVD)来帮助确定阶数。我们还可以指定一个阶数范围,然后选择使得信息准则最小的模型。 步骤: 1. 将数据组合成iddata对象。 2. 指定阶数范围(例如1:10)。 3. 使用`n4sid`命令估计不同阶数的模型,并返回每个模型的信息准则(如AIC)。 4. 选择AIC最小的模型对应的阶数。 注意:我们也可以使用`ssest`,但`n4sid`对于阶数选择更直观。 代码步骤详解: - 创建iddata对象:将输入u1,u2合并为一个矩阵,输出为y,指定采样时间(如果未知,可以设为1)。 - 使用`n4sidOptions`设置估计选项,例如选择使用AIC准则作为终止条件之一,但更常见的是我们指定阶数范围,然后计算每个阶数的模型,再比较AIC。 - 使用循环或直接使用`n4sid`的阶数范围参数来生成模型数组。 - 计算每个模型的AIC(或BIC)值,并找到最小值对应的阶数。 但是,MATLAB的`n4sid`函数在指定阶数范围时,会返回一个模型数组,然后我们可以比较这些模型。 另外,我们还可以通过绘制奇异值(在子空间方法中)来帮助判断阶数,但这里我们主要使用信息准则。 下面给出完整代码及详细注释。 </think> 以下是使用MATLAB确定控制系统黑箱模型状态空间方程最佳阶数的完整代码及详细注释。代码基于输入信号$u_1, u_2$和输出信号$y$,通过信息准则和模型验证选择最优阶数。 ```matlab %% 控制系统黑箱模型最佳阶数确定 % 作者: MATLAB助手 % 日期: 2023-10-15 % 输入: u1, u2 - 输入信号向量, y - 输出信号向量 % 输出: 最佳系统阶数n_optimal及对应模型 %% 步骤1: 数据准备与预处理 clc; clear; close all; % 加载示例数据(实际使用时替换为您的数据) load('systemData.mat'); % 假设数据文件包含u1, u2, y变量 % 检查数据维度 if size(u1, 2) > 1, u1 = u1'; end % 转为列向量 if size(u2, 2) > 1, u2 = u2'; end if size(y, 2) > 1, y = y'; end % 合并输入信号 (MISO系统) U = [u1, u2]; % 输入矩阵[N×2] % 创建iddata对象 (系统辨识工具箱) Ts = 0.1; % 采样时间(根据实际情况修改) data = iddata(y, U, Ts, 'Name', 'ControlSystem'); data.InputName = {'u1', 'u2'}; data.OutputName = {'y'}; % 数据拆分: 70%训练, 30%验证 N = size(y, 1); splitIdx = floor(0.7*N); trainData = data(1:splitIdx); valData = data(splitIdx+1:end); %% 步骤2: 阶数范围设置与初始化 n_min = 1; % 最小阶数 n_max = 10; % 最大阶数(根据系统复杂度调整) n_range = n_min:n_max; num_orders = length(n_range); % 初始化存储矩阵 models = cell(num_orders, 1); % 存储模型 AIC_values = zeros(num_orders, 1); % AIC值 BIC_values = zeros(num_orders, 1); % BIC值 MSE_values = zeros(num_orders, 1); % 均方误差 %% 步骤3: 模型估计与评估 (核心循环) fprintf('开始阶数搜索 (%d-%d)...\n', n_min, n_max); for i = 1:num_orders n = n_range(i); fprintf('拟合阶数 n=%d...', n); % 使用N4SID算法估计状态空间模型 options = n4sidOptions('Focus', 'prediction', 'N4Weight', 'auto'); sys = n4sid(trainData, n, 'Form', 'canonical', 'DisturbanceModel', 'none', options); % 计算验证集预测误差 y_pred = predict(sys, valData, 1); % 1步预测 residual = valData.y - y_pred.y; MSE = mean(residual.^2); % 计算信息准则 [~, AIC, BIC] = computeCriteria(sys, residual, n); % 保存结果 models{i} = sys; MSE_values(i) = MSE; AIC_values(i) = AIC; BIC_values(i) = BIC; fprintf('完成! MSE=%.4f, AIC=%.2f, BIC=%.2f\n', MSE, AIC, BIC); end %% 步骤4: 最佳阶数选择 (基于信息准则) % AIC准则: 选择最小值 [~, idx_AIC] = min(AIC_values); n_optimal_AIC = n_range(idx_AIC); % BIC准则: 选择最小值 [~, idx_BIC] = min(BIC_values); n_optimal_BIC = n_range(idx_BIC); % 综合决策 (优先BIC,因其惩罚项更强) if n_optimal_AIC == n_optimal_BIC n_optimal = n_optimal_AIC; fprintf('\n最佳阶数: n = %d (AIC和BIC一致)\n', n_optimal); else % 当结果不一致时,选择MSE变化率饱和点 MSE_ratio = abs(diff(MSE_values))./MSE_values(1:end-1); MSE_change_threshold = find(MSE_ratio < 0.05, 1); % MSE变化<5% n_optimal = n_range(min([idx_AIC, idx_BIC, MSE_change_threshold])); fprintf('\n最佳阶数: n = %d (AIC:%d, BIC:%d, MSE饱和点:%d)\n', ... n_optimal, n_optimal_AIC, n_optimal_BIC, n_range(MSE_change_threshold)); end %% 步骤5: 模型验证与可视化 % 绘制准则变化曲线 figure('Name', '阶数选择分析', 'Position', [100, 100, 1200, 500]) subplot(1,3,1) plot(n_range, MSE_values, 'bo-', 'LineWidth', 1.5) xlabel('系统阶数n'); ylabel('均方误差(MSE)'); title('MSE随阶数变化'); grid on; hold on; plot(n_optimal, MSE_values(n_optimal == n_range), 'ro', 'MarkerSize', 10); subplot(1,3,2) plot(n_range, AIC_values, 'go-', 'LineWidth', 1.5) xlabel('系统阶数n'); ylabel('AIC值'); title('AIC准则'); grid on; hold on; plot(n_optimal, AIC_values(n_optimal == n_range), 'ro', 'MarkerSize', 10); subplot(1,3,3) plot(n_range, BIC_values, 'mo-', 'LineWidth', 1.5) xlabel('系统阶数n'); ylabel('BIC值'); title('BIC准则'); grid on; hold on; plot(n_optimal, BIC_values(n_optimal == n_range), 'ro', 'MarkerSize', 10); % 残差相关性分析 (白噪声检验) best_model = models{n_optimal == n_range}; figure('Name', '残差分析') resid(best_model, trainData); % 残差自相关检验 title(sprintf('n=%d模型残差自相关函数', n_optimal)); % 模型响应验证 figure('Name', '模型验证') compare(valData, best_model); % 比较预测输出与实际输出 legend('真实输出', '模型预测', 'Location', 'best'); title(sprintf('n=%d模型验证集表现', n_optimal)); %% 步骤6: 输出最佳模型参数 fprintf('\n=== 最佳模型参数 (n=%d) ===\n', n_optimal); disp('A矩阵(状态转移矩阵):'); disp(best_model.A); disp('B矩阵(输入矩阵):'); disp(best_model.B); disp('C矩阵(输出矩阵):'); disp(best_model.C); disp('D矩阵(直接传输矩阵):'); disp(best_model.D); %% 辅助函数: 计算信息准则 function [FPE, AIC, BIC] = computeCriteria(sys, residuals, n) % 计算模型评价准则 % 输入: sys - 模型对象, residuals - 残差向量, n - 系统阶数 % 输出: FPE - 最终预测误差, AIC - Akaike信息准则, BIC - 贝叶斯信息准则 N = length(residuals); % 样本数量 V = det(cov(residuals)); % 残差协方差矩阵的行列式 d = size(sys, 'nu') * n + size(sys, 'ny') * n; % 参数总数 % 准则计算 FPE = V * (1 + d/N) / (1 - d/N); % 最终预测误差 AIC = log(V) + 2*d/N; % Akaike信息准则 BIC = log(V) + d*log(N)/N; % 贝叶斯信息准则 end ``` ### 代码说明 #### 1. **数据准备** - 输入信号合并为$N×2$矩阵$U = [u_1, u_2]$ - 创建`iddata`对象封装系统数据 - 数据分割为训练集(70%)和验证集(30%) #### 2. **阶数搜索设置** - 设置阶数范围$n_{\min}=1$到$n_{\max}=10$ - 初始化存储矩阵记录模型性能指标 #### 3. **模型估计循环** - 使用`n4sid`算法估计规范型状态空间模型 - 核心状态空间方程: $$ \begin{align*} \mathbf{x}(k+1) &= \mathbf{A}\mathbf{x}(k) + \mathbf{B}\mathbf{u}(k) \\ \mathbf{y}(k) &= \mathbf{C}\mathbf{x}(k) + \mathbf{D}\mathbf{u}(k) \end{align*} $$ - 计算验证集均方误差: $\text{MSE} = \frac{1}{N_v}\sum{(y-\hat{y})^2}$ - 计算信息准则(辅助函数`computeCriteria`): - $\text{AIC} = \ln(V) + \frac{2d}{N}$ - $\text{BIC} = \ln(V) + \frac{d\ln(N)}{N}$ #### 4. **最佳阶数选择** - 分别选择AIC和BIC最小的阶数 - 当结果不一致时,选择MSE变化率<5%的饱和点 - 决策公式: $\left|\frac{\Delta \text{MSE}}{\text{MSE}}\right| < 0.05$ #### 5. **模型验证** - 残差自相关函数检验(白噪声测试) - 模型预测输出与实际输出比较 - 参数矩阵$(\mathbf{A,B,C,D})$显示 ### 使用说明 1. **数据准备**:将实际输入输出数据赋值给`u1,u2,y` 2. **参数调整**: - 修改`Ts`为实际采样时间 - 调整`n_min`和`n_max`改变搜索范围 3. **结果解读**: - 选择准则曲线"拐点"对应的阶数 - 残差自相关函数应在置信区间内波动 - 验证集预测精度应满足工程需求 > **注意事项**: > 1. 需安装MATLAB系统辨识工具箱(System Identification Toolbox) > 2. 当数据噪声较大时,可增加`n_max`值 > 3. 对于复杂系统,建议结合奇异值分析法验证[^1]
评论 4
成就一亿技术人!
拼手气红包6.0元
还能输入1000个字符
 
红包 添加红包
表情包 插入表情
 条评论被折叠 查看
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

淘个代码_

不想刀我的可以选择爱我

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值