灰狼算法优化的LSSVR程序

使用灰狼优化算法(GWO)优化最小二乘支持向量回归(LSSVR)参数的MATLAB实现

classdef GWO_LSSVR
    % 灰狼算法优化的LSSVR模型
    % 使用灰狼优化算法寻找LSSVR的最佳参数(C, gamma)

    properties
        C;              % 正则化参数
        gamma;          % RBF核参数
        X_train;        % 训练数据特征
        y_train;        % 训练数据标签
        X_test;         % 测试数据特征
        y_test;         % 测试数据标签
        model;          % 训练好的LSSVR模型
        fitnessHistory; % 适应度历史记录
    end

    methods
        function obj = GWO_LSSVR(X_train, y_train, X_test, y_test)
            % 构造函数
            obj.X_train = X_train;
            obj.y_train = y_train;
            obj.X_test = X_test;
            obj.y_test = y_test;
            obj.fitnessHistory = [];
        end

        function [C_opt, gamma_opt] = train(obj, numWolves, maxIter)
            % 使用灰狼算法优化LSSVR参数
            % 参数设置
            dim = 2;                     % 优化参数维度 (C, gamma)
            lb = [0.1, 0.01];            % 参数下界 [C_min, gamma_min]
            ub = [1000, 100];            % 参数上界 [C_max, gamma_max]

            % 初始化灰狼种群
            wolves = initializeWolves(numWolves, dim, lb, ub);
            alpha_pos = zeros(1, dim);   % α狼位置
            alpha_score = inf;           % α狼适应度
            beta_pos = zeros(1, dim);    % β狼位置
            beta_score = inf;            % β狼适应度
            delta_pos = zeros(1, dim);   % δ狼位置
            delta_score = inf;           % δ狼适应度

            % 迭代优化
            for iter = 1:maxIter
                a = 2 - iter * (2 / maxIter);  % 线性递减参数a

                for i = 1:numWolves
                    % 边界检查
                    wolves(i, :) = boundConstraint(wolves(i, :), lb, ub);

                    % 计算适应度
                    fitness = obj.evaluateFitness(wolves(i, 1), wolves(i, 2));

                    % 更新α, β, δ狼
                    if fitness < alpha_score
                        alpha_score = fitness;
                        alpha_pos = wolves(i, :);
                    end

                    if fitness > alpha_score && fitness < beta_score
                        beta_score = fitness;
                        beta_pos = wolves(i, :);
                    end

                    if fitness > alpha_score && fitness > beta_score && fitness < delta_score
                        delta_score = fitness;
                        delta_pos = wolves(i, :);
                    end
                end

                % 更新每只狼的位置
                for i = 1:numWolves
                    for j = 1:dim
                        % 更新α狼影响
                        r1 = rand();
                        r2 = rand();
                        A1 = 2*a*r1 - a;
                        C1 = 2*r2;
                        D_alpha = abs(C1*alpha_pos(j) - wolves(i, j));
                        X1 = alpha_pos(j) - A1*D_alpha;

                        % 更新β狼影响
                        r1 = rand();
                        r2 = rand();
                        A2 = 2*a*r1 - a;
                        C2 = 2*r2;
                        D_beta = abs(C2*beta_pos(j) - wolves(i, j));
                        X2 = beta_pos(j) - A2*D_beta;

                        % 更新δ狼影响
                        r1 = rand();
                        r2 = rand();
                        A3 = 2*a*r1 - a;
                        C3 = 2*r2;
                        D_delta = abs(C3*delta_pos(j) - wolves(i, j));
                        X3 = delta_pos(j) - A3*D_delta;

                        % 位置更新
                        wolves(i, j) = (X1 + X2 + X3) / 3;
                    end
                end

                % 记录最佳适应度
                obj.fitnessHistory(iter) = alpha_score;

                % 显示迭代信息
                if mod(iter, 10) == 0
                    fprintf('Iteration %d: Best Fitness = %.4f, C = %.4f, gamma = %.4f\n', ...
                        iter, alpha_score, alpha_pos(1), alpha_pos(2));
                end
            end

            % 保存最优参数
            C_opt = alpha_pos(1);
            gamma_opt = alpha_pos(2);
            obj.C = C_opt;
            obj.gamma = gamma_opt;

            % 使用最优参数训练最终模型
            obj.model = obj.trainLSSVR(C_opt, gamma_opt);
        end

        function fitness = evaluateFitness(obj, C, gamma)
            % 评估适应度函数 (MSE)
            model = obj.trainLSSVR(C, gamma);
            y_pred = obj.predict(model, obj.X_test);
            mse = mean((obj.y_test - y_pred).^2);
            fitness = mse;  % 最小化MSE
        end

        function model = trainLSSVR(obj, C, gamma)
            % 训练LSSVR模型
            % 使用RBF核函数
            sigma = 1/sqrt(2*gamma);  % RBF核参数转换

            % 计算核矩阵
            n = size(obj.X_train, 1);
            K = zeros(n, n);
            for i = 1:n
                for j = 1:n
                    diff = obj.X_train(i, :) - obj.X_train(j, :);
                    K(i, j) = exp(-norm(diff)^2/(2*sigma^2));
                end
            end

            % 构建线性方程组
            Omega = [K + eye(n)/C, ones(n, 1); ones(n, 1)', 0];
            Y = [obj.y_train; 0];

            % 求解线性方程组
            sol = Omega \ Y;
            alpha = sol(1:n);
            b = sol(n+1);

            % 保存模型
            model.alpha = alpha;
            model.b = b;
            model.X_train = obj.X_train;
            model.gamma = gamma;
            model.C = C;
            model.sigma = sigma;
        end

        function y_pred = predict(obj, model, X)
            % 使用LSSVR模型进行预测
            n_train = size(model.X_train, 1);
            n_test = size(X, 1);
            y_pred = zeros(n_test, 1);

            for i = 1:n_test
                s = 0;
                for j = 1:n_train
                    diff = X(i, :) - model.X_train(j, :);
                    k_val = exp(-norm(diff)^2/(2*model.sigma^2));
                    s = s + model.alpha(j)*k_val;
                end
                y_pred(i) = s + model.b;
            end
        end

        function plotResults(obj)
            % 绘制优化过程和预测结果
            figure;

            % 适应度曲线
            subplot(2, 2, 1);
            plot(obj.fitnessHistory, 'LineWidth', 2);
            title('灰狼优化适应度曲线');
            xlabel('迭代次数');
            ylabel('均方误差(MSE)');
            grid on;

            % 预测结果与真实值对比
            y_pred = obj.predict(obj.model, obj.X_test);
            subplot(2, 2, 2);
            plot(obj.y_test, 'b-o', 'LineWidth', 1.5, 'DisplayName', '真实值');
            hold on;
            plot(y_pred, 'r-*', 'LineWidth', 1.5, 'DisplayName', '预测值');
            title('预测结果与真实值对比');
            xlabel('样本索引');
            ylabel('目标值');
            legend;
            grid on;

            % 误差分布
            subplot(2, 2, 3);
            errors = obj.y_test - y_pred;
            histogram(errors, 20);
            title('预测误差分布');
            xlabel('预测误差');
            ylabel('频数');
            grid on;

            % 参数优化结果
            subplot(2, 2, 4);
            text(0.1, 0.8, sprintf('最优参数 C = %.4f', obj.C), 'FontSize', 12);
            text(0.1, 0.6, sprintf('最优参数 gamma = %.4f', obj.gamma), 'FontSize', 12);
            text(0.1, 0.4, sprintf('最终MSE = %.4f', mean((obj.y_test - y_pred).^2)), 'FontSize', 12);
            text(0.1, 0.2, sprintf('R² = %.4f', 1 - sum((obj.y_test - y_pred).^2)/sum((obj.y_test - mean(obj.y_test)).^2)), 'FontSize', 12);
            axis off;
            title('优化结果总结');
        end

        function metrics = evaluateModel(obj)
            % 评估模型性能
            y_pred = obj.predict(obj.model, obj.X_test);

            % 计算性能指标
            mse = mean((obj.y_test - y_pred).^2);
            rmse = sqrt(mse);
            mae = mean(abs(obj.y_test - y_pred));
            r2 = 1 - sum((obj.y_test - y_pred).^2)/sum((obj.y_test - mean(obj.y_test)).^2);

            % 显示结果
            fprintf('\n===== 模型性能评估 =====\n');
            fprintf('均方误差(MSE): %.4f\n', mse);
            fprintf('均方根误差(RMSE): %.4f\n', rmse);
            fprintf('平均绝对误差(MAE): %.4f\n', mae);
            fprintf('决定系数(R²): %.4f\n', r2);

            % 返回指标
            metrics = struct('MSE', mse, 'RMSE', rmse, 'MAE', mae, 'R2', r2);
        end
    end
end

%% 辅助函数
function wolves = initializeWolves(numWolves, dim, lb, ub)
    % 初始化灰狼种群
    wolves = zeros(numWolves, dim);
    for i = 1:dim
        wolves(:, i) = lb(i) + (ub(i) - lb(i)) * rand(numWolves, 1);
    end
end

function x = boundConstraint(x, lb, ub)
    % 边界约束处理
    for i = 1:length(x)
        if x(i) < lb(i)
            x(i) = lb(i);
        elseif x(i) > ub(i)
            x(i) = ub(i);
        end
    end
end

%% 数据生成和演示脚本
function demoGWO_LSSVR()
    % 演示灰狼优化LSSVR
    clear; close all; clc;

    % 生成示例数据
    rng(42); % 设置随机种子
    nSamples = 200;
    nFeatures = 5;

    % 生成非线性关系数据
    X = rand(nSamples, nFeatures) * 10;
    y = sin(X(:,1)) + 0.5*cos(2*X(:,2)) + 0.3*X(:,3).^2 + 0.1*X(:,4).*X(:,5) + randn(nSamples, 1)*0.1;

    % 划分训练集和测试集
    splitRatio = 0.8;
    nTrain = round(splitRatio * nSamples);
    idx = randperm(nSamples);
    trainIdx = idx(1:nTrain);
    testIdx = idx(nTrain+1:end);

    X_train = X(trainIdx, :);
    y_train = y(trainIdx);
    X_test = X(testIdx, :);
    y_test = y(testIdx);

    % 创建并训练GWO-LSSVR模型
    gwo_lssvr = GWO_LSSVR(X_train, y_train, X_test, y_test);
    [C_opt, gamma_opt] = gwo_lssvr.train(20, 50); % 20只狼，50次迭代

    % 评估模型
    metrics = gwo_lssvr.evaluateModel();

    % 绘制结果
    gwo_lssvr.plotResults();

    % 与网格搜索比较
    compareWithGridSearch(X_train, y_train, X_test, y_test);
end

function compareWithGridSearch(X_train, y_train, X_test, y_test)
    % 与网格搜索调参比较
    fprintf('\n===== 网格搜索调参比较 =====\n');

    % 定义参数范围
    C_list = logspace(-1, 3, 10); % 0.1到1000
    gamma_list = logspace(-2, 2, 10); % 0.01到100

    best_mse = inf;
    best_C = 0;
    best_gamma = 0;

    % 网格搜索
    for i = 1:length(C_list)
        for j = 1:length(gamma_list)
            C = C_list(i);
            gamma = gamma_list(j);

            % 训练临时模型
            temp_model = trainLSSVR(X_train, y_train, C, gamma);
            y_pred = predictLSSVR(temp_model, X_test);
            mse = mean((y_test - y_pred).^2);

            if mse < best_mse
                best_mse = mse;
                best_C = C;
                best_gamma = gamma;
            end
        end
    end

    fprintf('网格搜索最佳参数: C = %.4f, gamma = %.4f\n', best_C, best_gamma);
    fprintf('网格搜索最佳MSE: %.4f\n', best_mse);

    % 使用最佳参数训练最终模型
    final_model = trainLSSVR(X_train, y_train, best_C, best_gamma);
    y_pred_grid = predictLSSVR(final_model, X_test);
    r2_grid = 1 - sum((y_test - y_pred_grid).^2)/sum((y_test - mean(y_test)).^2);
    fprintf('网格搜索R²: %.4f\n', r2_grid);
end

function model = trainLSSVR(X_train, y_train, C, gamma)
    % 训练LSSVR模型 (简化版)
    sigma = 1/sqrt(2*gamma);
    n = size(X_train, 1);
    K = zeros(n, n);
    for i = 1:n
        for j = 1:n
            diff = X_train(i, :) - X_train(j, :);
            K(i, j) = exp(-norm(diff)^2/(2*sigma^2));
        end
    end

    Omega = [K + eye(n)/C, ones(n, 1); ones(n, 1)', 0];
    Y = [y_train; 0];
    sol = Omega \ Y;

    model.alpha = sol(1:n);
    model.b = sol(n+1);
    model.X_train = X_train;
    model.sigma = sigma;
end

function y_pred = predictLSSVR(model, X)
    % 预测函数 (简化版)
    n_train = size(model.X_train, 1);
    n_test = size(X, 1);
    y_pred = zeros(n_test, 1);

    for i = 1:n_test
        s = 0;
        for j = 1:n_train
            diff = X(i, :) - model.X_train(j, :);
            k_val = exp(-norm(diff)^2/(2*model.sigma^2));
            s = s + model.alpha(j)*k_val;
        end
        y_pred(i) = s + model.b;
    end
end

程序功能说明

1. 核心类：GWO_LSSVR

• 属性：

  • C, gamma：LSSVR的正则化参数和核参数

  • X_train, y_train：训练数据

  • X_test, y_test：测试数据

  • model：训练好的LSSVR模型

  • fitnessHistory：优化过程中的适应度历史

• 方法：

  • train()：使用灰狼算法优化LSSVR参数

  • evaluateFitness()：评估给定参数的适应度(MSE)

  • trainLSSVR()：训练LSSVR模型

  • predict()：使用训练好的模型进行预测

  • plotResults()：可视化优化过程和结果

  • evaluateModel()：计算模型性能指标

2. 灰狼优化算法实现

• 初始化：在参数空间内随机生成狼群位置

• 迭代优化：

  • 更新参数a线性递减

  • 根据α、β、δ狼的位置更新其他狼的位置

  • 边界约束处理

• 适应度评估：使用测试集MSE作为适应度函数

3. LSSVR实现

• 核函数：径向基函数(RBF)

• 模型训练：求解线性方程组得到拉格朗日乘子和偏置

• 预测：使用核函数和学到的参数进行预测

4. 性能评估

• 计算多种回归指标：MSE、RMSE、MAE、R²

• 可视化预测结果与真实值对比

• 误差分布直方图

• 优化过程收敛曲线

5. 演示脚本

• 生成非线性回归数据集

• 划分训练集和测试集

• 训练GWO-LSSVR模型

• 与网格搜索调参方法进行性能比较

使用示例

% 运行演示程序
demoGWO_LSSVR();

% 自定义数据使用
% 1. 准备数据
load('your_dataset.mat'); % 加载数据
X = features; % 特征矩阵
y = targets;  % 目标变量

% 2. 划分数据集
splitRatio = 0.8;
nSamples = size(X, 1);
idx = randperm(nSamples);
nTrain = round(splitRatio * nSamples);
X_train = X(idx(1:nTrain), :);
y_train = y(idx(1:nTrain), :);
X_test = X(idx(nTrain+1:end), :);
y_test = y(idx(nTrain+1:end), :);

% 3. 创建并训练模型
gwo_lssvr = GWO_LSSVR(X_train, y_train, X_test, y_test);
[C_opt, gamma_opt] = gwo_lssvr.train(30, 100); % 30只狼，100次迭代

% 4. 评估模型
metrics = gwo_lssvr.evaluateModel();

% 5. 可视化结果
gwo_lssvr.plotResults();

参数调整建议

1. 灰狼算法参数：

   • numWolves：狼群数量（通常10-50）

   • maxIter：最大迭代次数（通常50-500）

   • 增加狼群数量和迭代次数可提高优化精度，但会增加计算时间

2. LSSVR参数范围：

   • lb和ub：参数搜索范围

   • 根据具体问题调整范围：

     • 对于简单问题：C∈[0.1, 100], gamma∈[0.01, 10]

     • 对于复杂问题：C∈[1, 1000], gamma∈[0.1, 100]

3. 核函数选择：

   • 当前实现使用RBF核

   • 可扩展其他核函数（多项式核、线性核等）

参考代码 使用灰狼算法优化的LSSVR程序 www.youwenfan.com/contentalh/100255.html

扩展功能建议

1. 多核学习支持

function model = trainMultiKernelLSSVR(obj, C, gamma1, gamma2, weight)
    % 训练多核LSSVR
    % 结合RBF核和多项式核
    n = size(obj.X_train, 1);
    K1 = zeros(n, n); % RBF核
    K2 = zeros(n, n); % 多项式核

    for i = 1:n
        for j = 1:n
            diff = obj.X_train(i, :) - obj.X_train(j, :);
            K1(i, j) = exp(-norm(diff)^2/(2*(1/sqrt(2*gamma1))^2));
            K2(i, j) = (1 + obj.X_train(i, :) * obj.X_train(j, :)') ^ gamma2;
        end
    end

    % 加权组合核矩阵
    K = weight * K1 + (1-weight) * K2;

    % 其余部分与单核相同...
end

2. 并行计算加速

function fitness = evaluateFitnessParallel(obj, wolves, C_list, gamma_list)
    % 并行评估适应度
    nWolves = size(wolves, 1);
    fitness = zeros(nWolves, 1);

    parfor i = 1:nWolves
        C = wolves(i, 1);
        gamma = wolves(i, 2);
        model = obj.trainLSSVR(C, gamma);
        y_pred = obj.predict(model, obj.X_test);
        fitness(i) = mean((obj.y_test - y_pred).^2);
    end
end

3. 混合优化策略

function [C_opt, gamma_opt] = hybridOptimization(obj, numWolves, maxIter)
    % 混合灰狼算法和局部搜索
    % 第一阶段：灰狼全局优化
    [C_gwo, gamma_gwo] = obj.train(numWolves, maxIter/2);

    % 第二阶段：局部搜索精调
    C_range = linspace(C_gwo*0.5, C_gwo*1.5, 10);
    gamma_range = linspace(gamma_gwo*0.5, gamma_gwo*1.5, 10);
    best_fitness = inf;

    for i = 1:length(C_range)
        for j = 1:length(gamma_range)
            fitness = obj.evaluateFitness(C_range(i), gamma_range(j));
            if fitness < best_fitness
                best_fitness = fitness;
                C_opt = C_range(i);
                gamma_opt = gamma_range(j);
            end
        end
    end
end

4. 时间序列预测应用

function forecast = timeSeriesForecast(obj, data, forecastHorizon)
    % 时间序列预测
    % data: 历史时间序列数据
    % forecastHorizon: 预测步长

    % 构造特征和目标
    X = [];
    y = [];
    for i = 1:(length(data)-obj.nFeatures)
        X(i, :) = data(i:i+obj.nFeatures-1);
        y(i) = data(i+obj.nFeatures);
    end

    % 训练模型
    obj.X_train = X;
    obj.y_train = y;
    obj.train(20, 50);

    % 递归预测
    lastValues = data(end-obj.nFeatures+1:end);
    forecasts = zeros(forecastHorizon, 1);

    for i = 1:forecastHorizon
        nextPred = obj.predict(obj.model, lastValues);
        forecasts(i) = nextPred;
        lastValues = [lastValues(2:end); nextPred];
    end

    forecast = forecasts;
end

实际应用建议

1. 数据预处理：

   • 标准化/归一化特征

   • 处理缺失值和异常值

   • 特征选择和降维（PCA等）

2. 参数调优：

   • 根据问题复杂度调整搜索范围

   • 使用敏感性分析确定关键参数

   • 考虑参数间的耦合关系

3. 模型集成：

   • 结合多个优化算法的优点

   • 集成多个LSSVR模型的预测

   • 使用bagging或boosting方法

4. 实时应用：

   • 增量学习更新模型

   • 滑动窗口处理新数据

   • 模型压缩减少计算开销