✅作者简介：热爱科研的Matlab仿真开发者，擅长毕业设计辅导、数学建模、数据处理、程序设计科研仿真。
🍎完整代码获取 定制创新 论文复现点击：Matlab科研工作室
👇 关注我领取海量matlab电子书和数学建模资料
🍊个人信条：做科研，博学之、审问之、慎思之、明辨之、笃行之，是为：博学慎思，明辨笃行。
🔥 内容介绍
一、引言
在科学与工程领域，混沌系统广泛存在，从气象学中的天气预测到电子电路中的信号处理，都能发现其身影。传统整数阶混沌系统已得到深入研究，但分数阶混沌系统由于其记忆性和遗传性，能更精确地描述许多复杂现象。基于数据驱动的分数阶混沌系统建模，通过对实际观测数据的分析，构建分数阶混沌模型，为理解和预测复杂动态系统提供了新的视角和方法。
二、分数阶混沌系统基础
（二）分数阶混沌系统特点
与整数阶混沌系统相比，分数阶混沌系统具有独特的动力学特性。由于分数阶导数的记忆性，分数阶混沌系统对初始条件的敏感性更强，其相空间轨迹更加复杂。例如，在一些分数阶混沌电路中，通过调整分数阶的阶数，可以观察到不同程度的混沌行为，且系统的吸引子形态也会发生变化。这种灵活性为模拟复杂的自然现象和工程系统提供了更强大的工具。
三、数据驱动建模的必要性与优势
（一）传统建模方法的局限
传统的混沌系统建模通常基于物理原理和数学推导，需要对系统的内部结构和参数有深入了解。然而，在许多实际场景中，系统的物理机制可能非常复杂，难以精确描述，或者系统参数难以直接测量。例如，在生物系统中，由于生物过程的高度复杂性和不确定性，基于物理原理的建模方法往往无法准确刻画系统的动态行为。
（二）数据驱动建模优势
- 无需精确物理模型
：数据驱动建模方法直接从观测数据出发，无需预先知道系统的物理模型。通过对大量数据的分析和挖掘，能够自动发现数据中的潜在规律，构建系统模型。这使得数据驱动建模适用于各种复杂系统，尤其是那些物理机制不明确的系统。
- 适应复杂动态变化
：实际系统往往处于动态变化中，数据驱动建模可以根据新的数据不断更新模型，适应系统的变化。例如，在电力系统中，负荷需求会随着时间、季节等因素变化，基于数据驱动的建模方法能够实时调整模型，准确预测电力负荷的变化。
四、基于数据驱动的分数阶混沌系统建模步骤
（一）数据采集
- 确定数据源
：根据研究对象，确定合适的数据采集方式和数据源。对于物理实验系统，可以通过传感器实时采集系统的输出数据，如在混沌电路实验中，采集电路的电压、电流等信号。对于一些实际应用场景，如气象数据，可以从气象站、卫星观测等渠道获取历史数据。
- 保证数据质量
：采集到的数据应具有足够的长度和精度，以反映系统的动态特性。同时，要对数据进行预处理，去除噪声和异常值。可以采用滤波算法（如卡尔曼滤波）对数据进行平滑处理，提高数据的质量。
（二）特征提取
- 混沌特征分析
：利用混沌理论中的方法，如相空间重构、李雅普诺夫指数计算等，分析采集数据的混沌特性。相空间重构通过延迟坐标法将一维时间序列数据映射到高维空间，恢复系统的动力学信息。李雅普诺夫指数则用于衡量系统对初始条件的敏感性，正的李雅普诺夫指数是混沌系统的重要特征之一。
- 分数阶特征挖掘
：从数据中挖掘与分数阶相关的特征。由于分数阶系统具有记忆性，数据的自相关性在分数阶特征提取中具有重要作用。可以通过分析数据的自相关函数在不同时间尺度上的衰减特性，提取与分数阶相关的信息。
（三）模型构建
- 选择建模方法
：常见的数据驱动建模方法包括神经网络、支持向量机回归等。神经网络具有强大的非线性拟合能力，能够逼近任意复杂的函数关系。支持向量机回归则在小样本数据建模中表现出色，能够有效避免过拟合问题。根据数据特点和建模需求，选择合适的建模方法。
- 确定模型结构
：以神经网络为例，需要确定网络的层数、每层的神经元个数等结构参数。可以通过交叉验证、网格搜索等方法，优化模型结构，提高模型的泛化能力。对于分数阶混沌系统建模，还需要考虑如何将分数阶微积分引入模型中，例如，通过在神经网络的隐藏层中加入分数阶微分算子，构建分数阶神经网络模型。
（四）模型验证与优化
- 模型验证
：将采集的数据划分为训练集、验证集和测试集。使用训练集对模型进行训练，利用验证集调整模型的超参数，最后在测试集上验证模型的性能。常用的验证指标包括均方误差（MSE）、平均绝对误差（MAE）等，这些指标用于衡量模型预测值与实际值之间的偏差。
- 模型优化
：根据验证结果，对模型进行优化。如果模型在测试集上的误差较大，可以尝试调整模型结构、增加数据量或采用更复杂的特征提取方法。同时，也可以结合其他优化算法，如遗传算法、粒子群优化算法等，对模型的参数进行全局优化，提高模型的性能。
⛳️ 运行结果
📣 部分代码
function [t, y] = fde12(alpha,fdefun,t0,tfinal,y0,h,param,mu,mu_tol)
%FDE12 Solves an initial value problem for a non-linear differential
% equation of fractional order (FDE). The code implements the
% predictor-corrector PECE method of Adams-Bashforth-Moulton type
% described in [1].
%
% [T,Y] = FDE12(ALPHA,FDEFUN,T0,TFINAL,Y0,h) integrates the initial value
% problem for the FDE, or the system of FDEs, of order ALPHA > 0
% D^ALPHA Y(t) = FDEFUN(T,Y(T))
% Y^(k)(T0) = Y0(:,k+1), k=0,...,m-1
% where m is the smallest integer grater than ALPHA and D^ALPHA is the
% fractional derivative according to the Caputo's definition. FDEFUN is a
% function handle corresponding to the vector field of the FDE and for a
% scalar T and a vector Y, FDEFUN(T,Y) must return a column vector. The
% set of initial conditions Y0 is a matrix with a number of rows equal to
% the size of the problem (hence equal to the number of rows of the
% output of FDEFUN) and a number of columns depending on ALPHA and given
% by m. The step-size H>0 is assumed constant throughout the integration.
%
% [T,Y] = FDE12(ALPHA,FDEFUN,T0,TFINAL,Y0,H,PARAM) solves as above with
% the additional set of parameters for the FDEFUN as FDEFUN(T,Y,PARAM).
%
% [T,Y] = FDE12(ALPHA,FDEFUN,T0,TFINAL,Y0,H,PARAM,MU) solves the FDE with
% the selected number MU of multiple corrector iterations. The following
% values for MU are admissible:
% MU = 0 : the corrector is not evaluated and the solution is provided
% just by the predictor method (the first order rectangular rule);
% MU > 0 : the corrector is evaluated by the selected number MU of times;
% the classical PECE method is obtained for MU=1;
% MU = Inf : the corrector is evaluated for a certain number of times
% until convergence of the iterations is reached (for convergence the
% difference between two consecutive iterates is tested).
% The defalut value for MU is 1
%
% [T,Y] = FDE12(ALPHA,FDEFUN,T0,TFINAL,Y0,H,PARAM,MU,MU_TOL) allows to
% specify the tolerance for testing convergence when MU = Inf. If not
% specified, the default value MU_TOL = 1.E-6 is used.
%
% FDE12 is an implementation of the predictor-corrector method of
% Adams-Bashforth-Moulton studied in [1]. Convergence and accuracy of
% the method are studied in [2]. The implementation with multiple
% corrector iterations has been proposed and discussed for multiterm FDEs
% in [3]. In this implementation the discrete convolutions are evaluated
% by means of the FFT algorithm described in [4] allowing to keep the
% computational cost proportional to N*log(N)^2 instead of N^2 as in the
% classical implementation; N is the number of time-point in which the
% solution is evaluated, i.e. N = (TFINAL-T)/H. The stability properties
% of the method implemented by FDE12 have been studied in [5].
%
% [1] K. Diethelm, A.D. Freed, The Frac PECE subroutine for the numerical
% solution of differential equations of fractional order, in: S. Heinzel,
% T. Plesser (Eds.), Forschung und Wissenschaftliches Rechnen 1998,
% Gessellschaft fur Wissenschaftliche Datenverarbeitung, Gottingen, 1999,
% pp. 57-71.
%
% [2] K. Diethelm, N.J. Ford, A.D. Freed, Detailed error analysis for a
% fractional Adams method, Numer. Algorithms 36 (1) (2004) 31-52.
%
% [3] K. Diethelm, Efficient solution of multi-term fractional
% differential equations using P(EC)mE methods, Computing 71 (2003), pp.
% 305-319.
%
% [4] E. Hairer, C. Lubich, M. Schlichte, Fast numerical solution of
% nonlinear Volterra convolution equations, SIAM J. Sci. Statist. Comput.
% 6 (3) (1985) 532-541.
%
% [5] R. Garrappa, On linear stability of predictor-corrector algorithms
% for fractional differential equations, Internat. J. Comput. Math. 87
% (10) (2010) 2281-2290.
%
% Copyright (c) 2011-2012, Roberto Garrappa, University of Bari, Italy
% garrappa at dm dot uniba dot it
% Revision: 1.2 - Date: July, 6 2012
% Check inputs
if nargin < 9
mu_tol = 1.0e-6 ;
if nargin < 8
mu = 1 ;
if nargin < 7
param = [] ;
end
end
end
% Check order of the FDE
if alpha <= 0
error('MATLAB:fde12:NegativeOrder',...
['The order ALPHA of the FDE must be positive. The value ' ...
'ALPHA = %f can not be accepted. See FDE12.'], alpha);
end
% Check the step--size of the method
if h <= 0
error('MATLAB:fde12:NegativeStepSize',...
['The step-size H for the method must be positive. The value ' ...
'H = %e can not be accepted. See FDE12.'], h);
end
% Structure for storing initial conditions
ic.t0 = t0 ;
ic.y0 = y0 ;
ic.m_alpha = ceil(alpha) ;
ic.m_alpha_factorial = factorial(0:ic.m_alpha-1) ;
% Structure for storing information on the problem
Probl.ic = ic ;
Probl.fdefun = fdefun ;
Probl.problem_size = size(y0,1) ;
Probl.param = param ;
% Check number of initial conditions
if size(y0,2) < ic.m_alpha
error('MATLAB:fde12:NotEnoughInputs', ...
['A not sufficient number of assigned initial conditions. ' ...
'Order ALPHA = %f requires %d initial conditions. See FDE12.'], ...
alpha,ic.m_alpha);
end
% Check compatibility size of the problem with size of the vector field
f_temp = f_vectorfield(t0,y0(:,1),Probl) ;
if Probl.problem_size ~= size(f_temp,1)
error('MATLAB:fde12:SizeNotCompatible', ...
['Size %d of the problem as obtained from initial conditions ' ...
'(i.e. the number of rows of Y0) not compatible with the ' ...
'size %d of the output of the vector field FDEFUN. ' ...
'See FDE12.'], Probl.problem_size,size(f_temp,1));
end
% Number of points in which to evaluate weights and solution
r = 16 ;
N = ceil((tfinal-t0)/h) ;
Nr = ceil((N+1)/r)*r ;
Q = ceil(log2(Nr/r)) - 1 ;
NNr = 2^(Q+1)*r ;
% Preallocation of some variables
y = zeros(Probl.problem_size,N+1) ;
fy = zeros(Probl.problem_size,N+1) ;
zn_pred = zeros(Probl.problem_size,NNr+1) ;
if mu > 0
zn_corr = zeros(Probl.problem_size,NNr+1) ;
else
zn_corr = 0 ;
end
% Evaluation of coefficients of the PECE method
nvett = 0 : NNr+1 ;
nalpha = nvett.^alpha ;
nalpha1 = nalpha.*nvett ;
PC.bn = nalpha(2:end) - nalpha(1:end-1) ;
PC.an = [ 1 , (nalpha1(1:end-2) - 2*nalpha1(2:end-1) + nalpha1(3:end)) ] ;
PC.a0 = [ 0 , nalpha1(1:end-2)-nalpha(2:end-1).*(nvett(2:end-1)-alpha-1)] ;
PC.halpha1 = h^alpha/gamma(alpha+1) ;
PC.halpha2 = h^alpha/gamma(alpha+2) ;
PC.mu = mu ; PC.mu_tol = mu_tol ;
% Initializing solution and proces of computation
t = t0 + (0 : N)*h ;
y(:,1) = y0(:,1) ;
fy(:,1) = f_temp ;
[y, fy] = Triangolo(1, r-1, t, y, fy, zn_pred, zn_corr, N, PC, Probl ) ;
% Main process of computation by means of the FFT algorithm
ff = [0 2 ] ; nx0 = 0 ; ny0 = 0 ;
for q = 0 : Q
L = 2^q ;
[y, fy] = DisegnaBlocchi(L, ff, r, Nr, nx0+L*r, ny0, t, y, fy, ...
zn_pred, zn_corr, N, PC, Probl ) ;
ff = [ff ff] ; ff(end) = 4*L ;
end
% Evaluation solution in TFINAL when TFINAL is not in the mesh
if tfinal < t(N+1)
c = (tfinal - t(N))/h ;
t(N+1) = tfinal ;
y(:,N+1) = (1-c)*y(:,N) + c*y(:,N+1) ;
end
t = t(1:N+1) ; y = y(:,1:N+1) ;
end
% =========================================================================
% =========================================================================
% r : dimension of the basic square or triangle
% L : factor of resizing of the squares
function [y, fy] = DisegnaBlocchi(L, ff, r, Nr, nx0, ny0, t, y, fy, ...
zn_pred, zn_corr, N , PC, Probl)
nxi = nx0 ; nxf = nx0 + L*r - 1 ;
nyi = ny0 ; nyf = ny0 + L*r - 1 ;
is = 1 ;
s_nxi(is) = nxi ; s_nxf(is) = nxf ; s_nyi(is) = nyi ; s_nyf(is) = nyf ;
i_triangolo = 0 ; stop = 0 ;
while ~stop
stop = nxi+r-1 == nx0+L*r-1 | (nxi+r-1>=Nr-1) ; % Ci si ferma quando il triangolo attuale finisce alla fine del quadrato
[zn_pred, zn_corr] = Quadrato(nxi, nxf, nyi, nyf, fy, zn_pred, zn_corr, PC, Probl) ;
[y, fy] = Triangolo(nxi, nxi+r-1, t, y, fy, zn_pred, zn_corr, N, PC, Probl) ;
i_triangolo = i_triangolo + 1 ;
if ~stop
if nxi+r-1 == nxf % Il triangolo finisce dove finisce il quadrato -> si scende di livello
i_Delta = ff(i_triangolo) ;
Delta = i_Delta*r ;
nxi = s_nxf(is)+1 ; nxf = s_nxf(is) + Delta ;
nyi = s_nxf(is) - Delta +1; nyf = s_nxf(is) ;
s_nxi(is) = nxi ; s_nxf(is) = nxf ; s_nyi(is) = nyi ; s_nyf(is) = nyf ;
else % Il triangolo finisce prima del quadrato -> si fa un quadrato accanto
nxi = nxi + r ; nxf = nxi + r - 1 ; nyi = nyf + 1 ; nyf = nyf + r ;
is = is + 1 ;
s_nxi(is) = nxi ; s_nxf(is) = nxf ; s_nyi(is) = nyi ; s_nyf(is) = nyf ;
end
end
end
end
% =========================================================================
% =========================================================================
function [zn_pred, zn_corr] = Quadrato(nxi, nxf, nyi, nyf, fy, zn_pred, zn_corr, PC, Probl)
coef_beg = nxi-nyf ; coef_end = nxf-nyi+1 ;
funz_beg = nyi+1 ; funz_end = nyf+1 ;
% Evaluation convolution segment for the predictor
vett_coef = PC.bn(coef_beg:coef_end) ;
vett_funz = [fy(:,funz_beg:funz_end) , zeros(Probl.problem_size,funz_end-funz_beg+1) ] ;
zzn_pred = real(FastConv(vett_coef,vett_funz)) ;
zn_pred(:,nxi+1:nxf+1) = zn_pred(:,nxi+1:nxf+1) + zzn_pred(:,nxf-nyf+1-1:end-1) ;
% Evaluation convolution segment for the corrector
if PC.mu > 0
vett_coef = PC.an(coef_beg:coef_end) ;
if nyi == 0 % Evaluation of the lowest square
vett_funz = [zeros(Probl.problem_size,1) , fy(:,funz_beg+1:funz_end) , zeros(Probl.problem_size,funz_end-funz_beg+1) ] ;
else % Evaluation of any square but not the lowest
vett_funz = [ fy(:,funz_beg:funz_end) , zeros(Probl.problem_size,funz_end-funz_beg+1) ] ;
end
zzn_corr = real(FastConv(vett_coef,vett_funz)) ;
zn_corr(:,nxi+1:nxf+1) = zn_corr(:,nxi+1:nxf+1) + zzn_corr(:,nxf-nyf+1:end) ;
else
zn_corr = 0 ;
end
end
% =========================================================================
% =========================================================================
function [y, fy] = Triangolo(nxi, nxf, t, y, fy, zn_pred, zn_corr, N, PC, Probl)
for n = nxi : min(N,nxf)
% Evaluation of the predictor
Phi = zeros(Probl.problem_size,1) ;
if nxi == 1 % Case of the first triangle
j_beg = 0 ;
else % Case of any triangle but not the first
j_beg = nxi ;
end
for j = j_beg : n-1
Phi = Phi + PC.bn(n-j)*fy(:,j+1) ;
end
St = StartingTerm(t(n+1),Probl.ic) ;
y_pred = St + PC.halpha1*(zn_pred(:,n+1) + Phi) ;
f_pred = f_vectorfield(t(n+1),y_pred,Probl) ;
if PC.mu == 0
y(:,n+1) = y_pred ;
fy(:,n+1) = f_pred ;
else
j_beg = nxi ;
Phi = zeros(Probl.problem_size,1) ;
for j = j_beg : n-1
Phi = Phi + PC.an(n-j+1)*fy(:,j+1) ;
end
Phi_n = St + PC.halpha2*(PC.a0(n+1)*fy(:,1) + zn_corr(:,n+1) + Phi) ;
yn0 = y_pred ; fn0 = f_pred ;
stop = 0 ; mu_it = 0 ;
while ~stop
yn1 = Phi_n + PC.halpha2*fn0 ;
mu_it = mu_it + 1 ;
if PC.mu == Inf
stop = norm(yn1-yn0,inf) < PC.mu_tol ;
if mu_it > 100 && ~stop
warning('MATLAB:fde12:NonConvegence',...
strcat('It has been requested to run corrector iterations until convergence but ', ...
'the process does not converge to the tolerance %e in 100 iteration'),PC.mu_tol) ;
stop = 1 ;
end
else
stop = mu_it == PC.mu ;
end
fn1 = f_vectorfield(t(n+1),yn1,Probl) ;
yn0 = yn1 ; fn0 = fn1 ;
end
y(:,n+1) = yn1 ;
fy(:,n+1) = fn1 ;
end
end
end
% =========================================================================
% =========================================================================
function z = FastConv(x, y)
r = length(x) ; problem_size = size(y,1) ;
z = zeros(problem_size,r) ;
X = fft(x,r) ;
for i = 1 : problem_size
Y = fft(y(i,:),r) ;
Z = X.*Y ;
z(i,:) = ifft(Z,r) ;
end
end
% =========================================================================
% =========================================================================
function f = f_vectorfield(t,y,Probl)
if isempty(Probl.param)
f = feval(Probl.fdefun,t,y) ;
else
f = feval(Probl.fdefun,t,y,Probl.param) ;
end
end
% =========================================================================
% =========================================================================
function ys = StartingTerm(t,ic)
ys = zeros(size(ic.y0,1),1) ;
for k = 1 : ic.m_alpha
ys = ys + (t-ic.t0)^(k-1)/ic.m_alpha_factorial(k)*ic.y0(:,k) ;
end
end
