使用 LSSVM 的 Matlab 演示求解反常微分方程问题（Matlab代码实现）

       目录

💥1 概述

📚2 运行结果

🎉3 参考文献

👨‍💻4 Matlab代码

---

💥1 概述

LSSVM的特性

  1) 同样是对原始对偶问题进行求解，但是通过求解一个线性方程组（优化目标中的线性约束导致的）来代替SVM中的QP问题（简化求解过程），对于高维输入空间中的分类以及回归任务同样适用；

  2) 实质上是求解线性矩阵方程的过程，与高斯过程（Gaussian processes），正则化网络（regularization networks）和费雪判别分析（Fisher discriminant analysis）的核版本相结合；

  3) 使用了稀疏近似（用来克服使用该算法时的弊端）与稳健回归（稳健统计）；

  4) 使用了贝叶斯推断（Bayesian inference）；

  5) 可以拓展到非监督学习中：核主成分分析（kernel PCA）或密度聚类；

  6) 可以拓展到递归神经网络中。

📚2 运行结果

主函数部分代码：

% dot(x1)  = a * (x_2 -x_1)

% dot(x_2) = x_1 * (b- x_3) - x_2

% dot(x_3) = x_1 * x_2 -c* x_3

% 0 <=  t  < = t_f

% Initial Condition

% x_1(0) = -9.42, x_2(0)= -9.34, x_3(0)=28.3

% Theta=[a, b, c] = [10, 28, 8/3]

%% ============================================================================

clear all; close all; clc

t0=0;

tf=10;

sampling_time=0.05;

t=(t0:sampling_time:tf)';

initial=[-9.42 -9.34 28.3]; % initial values of the ODE used for generating simulated data

ExactTheta=[10; 28 ; 8/3];  % The exact parameters of the lorenz system used for generating simulated data

cprintf( [1 0.1 0],'**** Excat parameters of the Lorenz system ***** \n\n');

fprintf('True theta_1= %f \n', ExactTheta(1));

fprintf('True theta_2= %f \n', ExactTheta(2));

fprintf('True theta_3= %f \n\n', ExactTheta(3));

fprintf( '************************************* \n\n');

%%  ========= Generating the simulation data ======================

options = odeset('RelTol',1e-5,'AbsTol',[1e-5 1e-5 1e-5]);

sol = ode45(@ridg,[t0 tf],initial,options,ExactTheta);

Y=deval(sol,t);

Y=Y';

noise_level=0.01; % 0.03, 0.05, 0.07, 0.1

noise=noise_level*randn(size(t,1),1);

y1=Y(:,1)+noise;

y2=Y(:,2)+noise;

y3=Y(:,3)+noise;

%% Estimating the parameters of the ODE system:

num_realization =3;

K_fold=3;

num_grid_gam=10;

num_grid_sig=10;

gamma_range = logspace(0,6,num_grid_gam);

sigma_range = logspace(-3,1,num_grid_sig);

ER1=[];

ER2=[];

ER3=[];

Par1=zeros(num_realization,1);

Par2=zeros(num_realization,1);

Par3=zeros(num_realization,1);

BB1=zeros(num_grid_gam,num_grid_sig);

BB2=zeros(num_grid_gam,num_grid_sig);

BB3=zeros(num_grid_gam,num_grid_sig);

for itr=1:num_realization

    

    

    cprintf( [1 0.1 0],'**** iteration =%d\n',itr);

    n=size(t,1);

    ind=crossvalind('Kfold', n, K_fold);

    

    for gamma_idx=1:size(gamma_range,2)

        gamma = gamma_range(gamma_idx);

        

        for sig_idx=1:size(sigma_range,2)

            sig = sigma_range(sig_idx);

            

            test_errors_1=zeros(K_fold,1);

            test_errors_2=zeros(K_fold,1);

            test_errors_3=zeros(K_fold,1);

            

            

            for i=1:K_fold

                Xte=t(ind==i,1);

                Yte_1=y1(ind==i,1);

                Yte_2=y2(ind==i,1);

                Yte_3=y3(ind==i,1);

                Xtr=t(ind~=i,1);

                Ytr_1=y1(ind~=i,1);

                Ytr_2=y2(ind~=i,1);

                Ytr_3=y3(ind~=i,1);

                               

                K=KernelMatrix(Xtr,'RBF_kernel', sig);

                m=size(K,1);

                A= [K + (1/gamma) * eye(m), ones(m,1);...

                    ones(m,1)' ,0];

                B1= [Ytr_1;0];

                B2= [Ytr_2;0];

                B3= [Ytr_3;0];

                result1=A\B1;

                result2=A\B2;

                result3=A\B3;

                alpha1=result1(1:m);

                b1=result1(end);

                alpha2=result2(1:m);

                b2=result2(end);

                alpha3=result3(1:m);

                b3=result3(end);

                yhattr1 = K * alpha1 + b1;

                yhattr2 = K * alpha2 + b2;

                yhattr3 = K * alpha3 + b3;

🎉3 参考文献

[1]姜星宇. 基于动态粒子群算法的DPSO-LSSVM模型在短期电力负荷预测中的应用研究[D].沈阳农业大学,2022.DOI:10.27327/d.cnki.gshnu.2022.000596.

👨‍💻4 Matlab代码