连续时间混沌系统MATLAB程序和SIMULINK模型

第6章连续时间混沌系统

本章讨论连续时间混沌系统的基本特点与分析方法，主要包括混沌数值仿真和硬件实验方法简介、混沌系数平衡点的计算、平衡点的分类与性质、相空间中的轨道、几类典型连续混沌系统的介绍、混沌机理的分析方法、用特征向量空间法寻找异宿轨道、Lorenz系统及混沌机理定性分析、Lorenz映射、Poincare截面、Chua系统及其混沌机理定性分析、时间序列与相空间重构等内容。

6.1 混沌数值仿真和硬件实验方法简介

混沌的数值仿真主要包括MA TLAB编程、SIMULINK模块构建、EWB仿真以及其他一些相关的软件仿真或数值计算等方法，从而获取混沌吸引子的相图、时域波形图、李氏指数、分叉图和功率谱等。混沌的硬件实验主要包括模拟/数字电路设计与硬件实验、现场可编程门阵列器件（FPGA）、数字信号处理器（DSP）等硬件实现方法来产生混沌信号。本节仅对各种数值仿真方法作简单介绍。

1）混沌系统的MA TLAB数值仿真

该方法主要根据混沌系统的状态方程来编写MA TLAB程序。现举二例来说明这种编程方法。（1）已知Lorenz系统的状态方程为

dx/dt=-a(x-y)

dy/dt=bx-xz-y

dz/dt=-cz+xy

式中a=10,b=30,c=8/3。

MA TLAB仿真程序如下：

>> %**************************************************

Function dxdt=lorenz(t,x)

%除符号dxdt外，还可用其他编程者习惯的有意义的符号

A=10;

B=30;

C=8/3;

dxdt=zeros(3,1);

dxdt(1)=-A*(x(1)-x(2));

dxdt(2)=B*x(1)-x(1).*x(3)-x(2);

dxdt(3)=x(1)*x(2)-C*x(3);

%*************************************************

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

t0=[0 200];

x0=[0.02,0.01,0.03];

[t,x]=ode45('lorenz',t0,x0,options);

%**************************************************

n=length(t)

n1=round(n/2)

%n1=1;

%**************************************************

figure(1);

plot(t(n1:n,1),x(n1:n,1));

xlabel('t','fontsize',20,'fontname','times new roman','FontAngle','italic');

ylabel('x','fontsize',20,'fontname','times new roman','FontAngle','italic');

figure(2);

plot(x(n1:n,1),x(n1:n,3));

xlabel('x','fontsize',20,'fontname','times new roman','FontAngle','italic');

ylabel('z','fontsize',20,'fontname','times new roman','FontAngle','italic');

%*******************************************************************

根据上述MA TLAB程序，得Lorenz系统的时域波形图和混沌吸引子相图的数值仿真结果如图6-1所示。

图6-1 Lorenz系统的时域波形图和混沌吸引子相图的MATLAB数值仿真结果

（2）已知Chua系统的状态方程为

dx=a[y-f(x)]

dy=x-y+z

dz=-by

式中a=10,b=15,m0=-1/7,m1=2/7,f(x)=m1*x+0.5(m0-m1)[|x+1|-|x-1|]为三分段非线性函数

MA TLAB仿真如下：

function dxdt=chua(t,x)

m0=-1/7;m1=2/7;

a=10;

b=15;

%*******************************************

dxdt=zeros(3,1);

f=m1*x(1)+0.5*(m0-m1)*(abs(x(1)+1)-abs(x(1)-1));

dxdt(1)=a*(x(2)-f);

dxdt(2)=x(1)-x(2)+x(3);

dxdt(3)=-b*x(2);

%*****************************************************

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

t0=[0 5e+2];

x0=[0.01 0.02 0.03];

[t,x]=ode45('chua',t0,x0,options);

%****************************************************

n=length(t)

n1=round(n/2)

%****************************************************** figure(1);

plot(t(n1:n),x(n1:n,1));

xlabel('t','fontsize',20,'fontname','times new roman','FontAngle','italic'); ylabel('x','fontsize',20,'fontname','times new roman','FontAngle','italic'); figure(2);

plot(x(n1:n,1),x(n1:n,2));

xlabel('x','fontsize',20,'fontname','times new roman','FontAngle','italic'); ylabel('y','fontsize',20,'fontname','times new roman','FontAngle','italic'); %********************************************************

根据上述MA TLAB 程序，得Chua 系统的时域波形图和混沌吸引子相图的数值仿真结果如图6-2所示。

图6-2 Chua 系统的时域波形图和混沌吸引子相图的MATLAB 数值仿真结果

2）混沌系统的SIMULINK 仿真

该方法主要是根据混沌系统的状态方程，将其转换成积分方程，利用模块和积分算子画出SIMULINK 的模块化仿真图。为保证计算的精确度，又不使仿真时间过长，应

对仿真图中几个重要参数进行设置。第一个参数是仿真时间：第二个参数是相对误差，通常设为10

10-；第三个参数是绝对误差，通常设为13

-，现举二例来说明这种编程

方法。

（1）已知Lorenz 系统的状态方程仍如（6-1）式，将其转换成积分方程：

dx/dt=-a(x-y) ?--=dt y x a x )]([

dy/dt=bx-xz-y ?--=

y xz bx

y ][ dz/dt=-cz+xy

?+-=

xy cz

z ][

注意，SIMULINK 仿真中的微分子算子为S ，积分算子为

1,故得SIMULINK 的仿真

如图6-3所示，设其文件名为“simulink_lorenz ”,再利用文件名为“y_simulink_lorenz ”的程序运行“simulink_lorenz ”。程序如下：

[t,x]=sim('simulink_lorenz',200);

n=length(t)

n1=round(n/2)

%*********************************************

figure(1);

plot(t(n1:n,1),x(n1:n,1));

xlabel('t','fontsize',20,'fontname','times new roman','FontAngle','italic');

ylabel('x','fontsize',20,'fontname','times new roman','FontAngle’,'italic');

figure(2);

plot(x(n1:n,1),x(n1:n,3));

xlabel('x','fontsize',20,'fontname','times new roman','FontAngle','italic');

ylabel('z','fontsize',20,'fontname','times new roman','FontAngle','italic');

%********************************************************

运行结果仍如图6-1所示。

图6-3 Lorenz系统的SIMULINK仿真图

（2）已知Chua系统的状态方程仍如（6-2）式，得SIMULINK的仿真图如图6-4所示。设其文件名为“simulink_chua”,再利用文件名为“y_simulink_chua”的程序运行“simulink_chua”。程序如下：

%*******************************************

%global a;

a=10;

[t,x]=sim('simulink_chua',500);

n=length(t)

n1=round(n/2)

%********************************************

figure(1);

plot(t(n1:n,1),x(n1:n,1));

xlabel('t','fontsize',20,'fontname','times new roman','FontAngel','italic');

ylabel('x','fontsize',20,'fontname','times new roman','FontAngel','italic');

figure(2);

plot(x(n1:n,1),x(n1:n,2));

xlabel('x','fontsize',20,'fontname','times new roman','FontAngel','italic'); ylabel('y','fontsize',20,'fontname','times new roman','FontAngel','italic');

运行结果仍如图6-2所示。

（a ）SIMULINK 主框图

（b ）三分段线性函数f 的SIMULINK 子框图

图6-4 Chua 系统的SIMULINK 仿真图

3）连续混沌系统离散化的MA TLAB 数值仿真

当用DSP 和FPGA 等现代数字器件来产生混沌信号时，首先需要将连续混沌系统作离散化处理。离散化和数字化处理方法主要有三种，利用这些离散化的方法，可将状态方程变成差分方程，这些方法将在后续章节中详细介绍。这里采用了一种较简单的Euler 算法。

（1）已知Lorenz 系统的状态方程仍如（6-1）式，根据Euler 算法

n x n x t

n x n x dt

dx )

()1()

()1(-+??→

??-+=

简记为

得对应的差分方程为

)

()()()

()1()()()()()

()1()]

()([)

()1(n y n x n cz T

n z n z n y n z n x n bx T

n y n y n y n x a T

n x n x +-=-+--=-+--=-+

将其整理成标准的差分方程形式：

)

()()()1()1()()1()()()()1()

()()1()1(n y n Tx n z Tc n z n y T n z n Tx n Tbx n y n Tay n x Ta n x +-=+-+-=++-=+

根据（6-3）式，得MA TLAB 仿真程序如下：

clear all; hold off;

%******************************************

x(1)=0.02;y(1)=0.01;z(1)=0.03; a=10;b=30;c=8/3; T=5e-3;%可调取样时间 N=5e4;%可调迭代次数 for i=1:N i t(i)=i;

x(i+1)=(1-a*T)*x(i)+a*T*y(i);

y(i+1)=T*b*x(i)-T*x(i)*z(i) +(1-T)*y(i); z(i+1)=(1-T*c)*z(i)+T*x(i)*y(i); end

%************************************ n1=N/2;n2=N; figure(1)

plot(t(n1:n2),x(n1:n2),’b-’);

xlabel('t_n','fontsize',20,'fontname','times new roman','FontAngle','italic'); ylabel('x_n','fontsize',20,'fontname','times new roman','FontAngle','italic'); figure(2);

plot(x(n1:n2),x(n1:n2),’b-’);

xlabel('x_n','fontsize',20,'fontname','times new roman','FontAngle','italic'); ylabel('z_n','fontsize',20,'fontname','times new roman','FontAngle','italic'); %************************************

运行结果仍如图6-1所示。

(2)已知Chua 系统的状态方程仍如（6-2）式，根据Euler 算法，得对应的差分方程为

)

()

()1()()()()

()1())]

(()([)

()1(n by T

n z n z n z n y n x T

n y n y n y f n y a T

n x n x -=-++-=-+-=-+

将其整理成标准的差分方程形式：

)

()()1()()()1()()1())

(()()()1(n z n Tb n z n Tz n y T n Tx n y n x Taf n Tay n x n x +-=++-+=+-+=+

式中a=10,b=15,f(x(n))=m1x(n)+0.5(m0-m1)[|x(n)+1|-|x(n)-1|]为离散化后的三分段线性函数，其中参数m0=-1/7,m1=2/7.根据（6-4）式，得MA TLAB 仿真程序如下：

clear all ; hold off ;

%****************************************** x(1)=0.02;y(1)=0.01;z(1)=0.03; m0=-1/7; m1=2/7; a=10;b=15;

T=1e-2;%?éμ÷è??ùê±?? N=5e4;%?éμ÷μü′ú′?êy for i=1:N i t(i)=i;

f(i)=m1*x(i)+0.5*(m0-m1).*(abs(x(i)+1)-abs(x(i)-1)); x(i+1)= x(i)+a*T*y(i)-T*a*f(i); y(i+1)=T *x(i)+(1-T)*y(i)+T*z(i); z(i+1)= -T*b*y(i)+z(i); end

%************************************ n1=N/2;n2=N; figure(1)

plot(t(n1:n2),x(n1:n2),'b-');

xlabel('t_n','fontsize',20,'fontname','times new roman','FontAngle','italic');

ylabel('x_n','fontsize',20,'fontname','times new roman','FontAngle','italic'); figure(2);

plot(x(n1:n2),y(n1:n2),'b-');

xlabel('x_n','fontsize',20,'fontname','times new roman','FontAngle','italic');

ylabel('y_n','fontsize',20,'fontname','times new roman','FontAngle','italic');

%************************************

运行结果仍如图6-2所示。

4）离散混沌系统的SIMULINK 仿真

仍以Lorenz 系统和Chua 系统为例来说明。 (1)Lorenz 系统离散化后的SIMULINK 仿真。为方便计，将（6-3）式表示为如下的标准形式：

)

()()()1()()()()()1()

()()1(n y n Gx n Fz n z n Ey n z n Dx n Cx n y n By n Ax n x +=+++=++=+

式中A=1-Ta,B=Ta,C=Tb,D=T,E=1-T,F=1-Tc,G=T,得离散化后Lorenz 系统的SIMULINK 仿真如图6-5所示。

注意，由于是离散系统，故应将连续系统中的算子1/S 换成离散系统中的算子1/Z 。设其文件为“simulink_discrete_lorenz ”的程序运行“simulink_discrete_lorenz ”。程序如下：

图6-5 离散化后Lorenz 系统的SIMULINK 仿真图

%********************************************* a=10;b=30;c=8/3;T=5e-3;%可调取样时间

%global A;global B;global D;global E;global F;global G; A=1-T*a;B=T*a;C=T*b;D=T;E=1-T*c;G=T;

%********************************************* [t,x]=sim(’simulink_discrete_lorenz ’,5e4); n=length(t) n1=round(n/2);

%********************************************* figure(1);

plot(t(n1:n,1),x(n1:n,1));

xlabel('t','fontsize',20,'fontname','times new roman','FontAngel','italic');

ylabel('x','fontsize',20,'fontname','times new roman','FontAngel','italic');

figure(2);

plot(x(n1:n,1),x(n1:n,2));

xlabel('x','fontsize',20,'fontname','times new roman','FontAngel','italic');

ylabel('z','fontsize',20,'fontname','times new roman','FontAngel','italic');

%*********************************************

运行结果仍如图6-1所示。

（2）Chua系统离散化后的SIMULINK仿真。为方便计，将（6-4）式表示为如下的标准形式：

)

(

)

(

) (

)

(

)

(

))

(

)

(

)

(

n Ey

x Cf

式中，A=1,B=Ta,C=Ta,D=T,E=1-T,F=T,G=Tb,H=1,T=1e-2,a=10,b=15,m0= -1/7,m1=2/7,得离散化后Lorenz系统的SIMULINK仿真图如6-6所示。

(a)SIMULINK主框图（b）三分段线性函数f(x(n))的SIMULINK子框图

图6-6 离散化后Chua系统的SIMULINK仿真图

设其文件名为“simulink_discrete_chua”,利用文件名为“y_simulink_discrete_chua”的程序运行“simulink_discrete_chua”。程序如下：

%*********************************************

a=10;b=15;T=1e-2;%可调取样时间

%global A;global B;global C;global D;

%global E;global F;global G;global H;

A=1;

B=T*a;

C=T*a;

D=T;

E=1-T;

G=T*b;

H=1;

%*********************************************

[t,x]=sim(’simulink_discrete_lorenz’,5e4);

n=length(t)

n1=round(n/2);

%*********************************************

figure(1);

plot(t(n1:n,1),x(n1:n,1));

xlabel('t','fontsize',20,'fontname','times new roman','FontAngel','italic');

ylabel('x','fontsize',20,'fontname','times new roman','FontAngel','italic');

figure(2);

plot(x(n1:n,1),x(n1:n,2));

xlabel('x','fontsize',20,'fontname','times new roman','FontAngel','italic');

ylabel('z','fontsize',20,'fontname','times new roman','FontAngel','italic');

%*********************************************

运行结果仍如图6-2所示。

以上我们分别以Lorenz系统和chua系统为典型实例，对四种数值仿真方法进行了详细的说明。这四种仿真方法分别是混沌系统的MA TLAB数值仿真、混沌系统的SIMULINK仿真、连续混沌系统离散化的MA TLAB数值仿真和离散混沌系统的SIMULINK仿真。有关用模拟电路、DSP、FPGA产生混沌的问题将在后续章节中给出，此处不再讨论。

蔡氏电路MATLAB混沌仿真

蔡氏电路的Matlab混沌仿真研究班级：姓名：学号：

摘要本文首先介绍非线性系统中的混沌现象，并从理论分析与仿真计算两个方面细致研究了非线性电路中典型混沌电路，即蔡氏电路反映出的非线性性质。通过改变蔡氏电路中元件的参数，进而产生多种类型混沌现象。最后利用软件对蔡氏电路的非线性微分方程组进行编程仿真，实现了双涡旋和单涡旋状态下的同步，并准确地观察到混沌吸引子的行为特征。关键词：混沌；蔡氏电路；MATLAB仿真 Abstract This paper introduce s the chaos phenomenon in nonlinear circuits. Chua’s circuit was a typical chaos circuit, thus theoretical analysis and simulation was made to research it. Many kinds of chaos phenomenon on would generate as long as one component parameter was altered in C hua’s circuit．On the platform of Matlab, mathematical model of Chua’s circuit was programmed and simulated to acquire the synchronization of dual and single cochlear volume. Meanwhile, behavioral characteristics of chaos attractor were observed. Key words：chaos phenomenon；Chua’s circuit；Simulation

用Matlab观察分岔与混沌现象

M a t l a b 实验报告实验目的：用Matlab 观察分岔与混沌现象。题目：Feigenbaum 曾对超越函数sin()y x λπ=（λ为非负实数）进行了分岔与混沌的研究，试利用迭代格式1sin()k k x x λπ+=，做出相应的Feigenbaum 图算法设计： 1、因为λ为非负实数，所以试将λ的范围限制在[0,3]，制图时x 的坐标限制在[0,3]，考虑到y 的值有正有负，所以把y 的坐标限制在 [-3,3]。 2、根据课本上给的例题，编写程序代码来绘图。程序代码： clear;clf; hold on axis([0,3,-3,3]); grid for a=0:0.005:3 x=[0.1]; for i=2:150 x(i)=a*sin(pi*x(i-1)); end pause(0.1) for i=101:150 plot(a,x(i),'k.'); end end 图像：结果分析：在λ取值在[0,0.3]区间内时，y 的值保持在0，然后开始上升，在λ取值在0.75附近时，开始分岔为两支。从整体上看，随着λ的值越来越大，所产生的迭代序列越来越复杂，可能会随机地落在区间（-3,3）的任一子区间内。并可能重复，这就是混沌的遍历性。进一步分析：由于λ的取值空间偏小，考虑扩大其取值范围

到[0,6]，再进一步观察图像。程序代码如下： clear;clf; hold on axis([0,6,-6,6]); grid for a=0:0.05:6 x=[0.1]; for i=2:150 x(i)=a*sin(pi*x(i-1)); end pause(0.1) for i=101:150 plot(a,x(i),'k.'); end end 图像：分析：由图像可见，随着取值范围的增大，图像呈现出周期性的特点。总结：1、当取值范围比较小，不足以发现图像规律时，可以考虑扩大变量的取值范围。 2、由于图像是由大量点构成的，所以在编程的时候注意循环语句的应用。

灰色预测模型的Matlab程序及检验程序(精)

灰色预测模型的Matlab 程序及检验程序 %灰色预测模型程序 clear syms a b; c=[a b]'; A=[46.2 32.6 26.7 23.0 20.0 18.9 17.5 16.3];% 原始序列 B=cumsum(A);%累加n=length(A); for i=1:(n-1) C(i)=(B(i)+B(i+1))/2; end %计算待定参数 D=A; D(1)=[]; D=D'; E=[-C; ones(1,n-1)]; c=inv(E*E')*E*D; c=c'; a=c(1); b=c(2); %预测往后预测5个数据 F=[];F(1)=A(1); for i=2:(n+5) F(i)=(A(1)-b/a)/exp(a*(i-1))+b/a; end G=[];G(1)=A(1); for i=2:(n+5) G(i)=F(i)-F(i-1); end t1=2002:2009; t2=2002:2014; G plot(t1,A,'o',t2,G) %灰色预测模型检验程序 function [ q,c,p ] = checkgm( x0,x1 ) %GM 检验函数 %x0 原始序列

%x1 预测序列 %·返回值 % q –- 相对误差 % c -- ·方差比 % p -- 小误差概率 e0=x0-x1; q=e0/x0; s1=var(x0); %qpa=mean(e0); s2=var(e0); c=s2/s1; len=length(e0); p=0; for i=1:len if(abs(e0(i)) < 0.6745*s1) p=p+1; end end p=p/len; end

(完整版)基于MATLAB的混沌序列图像加密程序

设计题目：基于MATLAB的混沌序列图像加密程序一.设计目的图像信息生动形象，它已成为人类表达信息的重要手段之一，网络上的图像数据很多是要求发送方和接受都要进行加密通信，信息的安全与保密显得尤为重要，因此我想运用异或运算将数据进行隐藏，连续使用同一数据对图像数据两次异或运算图像的数据不发生改变，利用这一特性对图像信息进行加密保护。熟练使用matlab运用matlab进行编程，使用matlab语言进行数据的隐藏加密，确保数字图像信息的安全，混沌序列具有容易生成，对初始条件和混沌参数敏感等特点，近年来在图像加密领域得到了广泛的应用。使用必要的算法将信息进行加解密，实现信息的保护。 .设计内容和要求使用混沌序列图像加密技术对图像进行处理使加密后的图像使用matlab将图像信息隐藏，实现信息加密。三.设计思路 1. 基于混沌的图像置乱加密算法本文提出的基于混沌的图像置乱加密算法示意图如图1所示加密算法如下：首先，数字图像B大小为MX N（ M是图像B的行像素数，N是图像B的列像素数），将A的第j行连接到j-1行后面（j=2，3, A，M，形成长度为MX N的序列C。其次，用Logistic混沌映射产生一个长度为的混沌序列｛k1，k2，A，kMX N｝，并构造等差序列D: ｛1，2，3, A，MX N-1，MX N｝。再次，将所

产生的混沌序列｛kl, k2. A, kMX N｝的M N个值由小到大排序，形成有序序列｛k1', k2'. A' kMX N' ｝，确定序列｛k1, k2, A, kMX N｝中的每个ki在有序序列｛k1', k2', A , kMX N' ｝中的编号，形成置换地址集合｛t1 , t2 , A, tM X N｝，其中ti为集合｛1 , 2, A, MX N｝中的一个；按置换地址集合｛t1 , t2 , A, tM X N｝对序列C进行置换，将其第i个像素置换至第ti列, i=1 , 2, A, MX N,得到C'。将等差序列D做相同置换，得到D'。最后，B'是一个MX N 的矩阵，B' (i ,j)=C ' ((i-1) X M+j)，其中i=1 , 2, A, M j=i=1 , 2, A, N,则B'就是加密后的图像文件。解密算法与加密算法相似，不同之处在于第3步中，以序列C'代替随机序列｛k1, k2, A, kMX N｝,即可实现图像的解密。 2. 用MATLAB勺实现基于混沌的图像置乱加密算法本文借助MATLAB^件平台，使用MATLAB!供的文本编辑器进行编程实现加密功能。根据前面加密的思路，把加密算法的编程分为三个主要模块：首先，构造一个与原图a等高等宽的矩阵b加在图像矩阵a后面形成复合矩阵c: b=zeros(m1, n1); ifm1>=n1 ifm1> n1 fore=1: n1 b=(e,e); end else fore=1: n1 end fore=1:( n1-m1) b((m1+e-1),e)=m1+e-1 end end c=zeros(m1*2, n1); c=zeros(m1*2,1); c=[b,a]; 然后，用Logitic映射产生混沌序列:

灰色预测MATLAB程序

作用：求累加数列、求a b的值、求预测方程、求残差 clc %清屏，以使结果独立显示 x=[ ]; format long; %设置计算精度 if length(x(:,1))==1 %对输入矩阵进行判断，如不是一维列矩阵，进行转置变换 x=x'; end n=length(x); %取输入数据的样本量 z=0; for i=1:n %计算累加值，并将值赋予矩阵be z=z+x(i,:); be(i,:)=z; end for i=2:n %对原始数列平行移位 y(i-1,:)=x(i,:); end

for i=1:n-1 %计算数据矩阵B的第一列数据 c(i,:)=*(be(i,:)+be(i+1,:)); end for j=1:n-1 %计算数据矩阵B的第二列数据 e(j,:)=1; end for i=1:n-1 %构造数据矩阵B B(i,1)=c(i,:); B(i,2)=e(i,:); end alpha=inv(B'*B)*B'*y; %计算参数矩阵即a b的值 for i=1:n+1 %计算数据估计值的累加数列，如改为n+1为n+m可预测后m-1个值 ago(i,:)=(x(1,:)-alpha(2,:)/alpha(1,:))*exp(-alpha(1,:)*(i-1))+alpha(2,:)/alpha(1,: );%显示输出预测值的累加数列 end var(1,:)=ago(1,:) %显示输出预测值 for i=1:n %如改n为n+m-1，可预测后m-1个值 var(i+1,:)=ago(i+1,:)-ago(i,:); %估计值的累加数列的还原，并计算出下一预测值end for i=1:n error(i,:)=x(i,:)-var(i,:); %计算残差 end c=std(error)/std(x); %调用统计工具箱的标准差函数计算后验差的比值c ago %显示输出预测值的累加数列 alpha %显示输出参数数列 var %显示输出预测值 error %显示输出误差 c %显示后验差的比值作用：数据处理判断是否可以用灰色预测、求级比、求累加数列、求a b的值、求预测方程 clc,clear x0=[ ]'; %注意这里为列向量 n=length(x0); lamda=x0(1:n-1)./x0(2:n) %计算级比 range=minmax(lamda') %计算级比的范围 x1=cumsum(x0) %累加运算 B=[*(x1(1:n-1)+x1(2:n)),ones(n-1,1)]; Y=x0(2:n); u=B\Y %拟合参数u(1)=a,u(2)=b x=dsolve('Dx+a*x=b','x(0)=x0'); %求微分方程的符号解

灰色预测matlab程序

Matlab的灰色预测程序： y=input('请输入数据'); n=length(y); yy=ones(n,1); yy(1)=y(1); for i=2:n yy(i)=yy(i-1)+y(i) end B=ones(n-1,2); for i=1:(n-1) B(i,1)=-(yy(i)+yy(i+1))/2; B(i,2)=1; end BT=B'; for j=1:(n-1) YN(j)=y(j+1); end YN=YN'; A=inv(BT*B)*BT*YN; a=A(1); u=A(2); t=u/a; t_test=input('输入需要预测的个数'); i=1:t_test+n; yys(i+1)=(y(1)-t).*exp(-a.*i)+t; yys(1)=y(1); for j=n+t_test:-1:2 ys(j)=yys(j)-yys(j-1); end x=1:n; xs=2:n+t_test; yn=ys(2:n+t_test); plot(x,y,'^r',xs,yn,'*-b'); det=0; for i=2:n det=det+abs(yn(i)-y(i)); end det=det/(n-1); disp(['百分绝对误差为：',num2str(det),'%']); disp(['预测值为：',num2str(ys(n+1:n+t_test))]); 请输入数据[29.8 30.11 41.05 70.12 77.79 77.79 104.82 65.22 82.7 100.79] 输入需要预测的个数4

百分绝对误差为：14.5128% 预测值为：110.5718 120.8171 132.0116 144.2434

灰色预测模型matlab程序精确版

灰色预测模型matlab程序 %下面程序是灰色模型GM(1,1)程序二次拟合和等维新陈代谢改进预测程序,mat lab6.5 ,使用本程序请注明，程序存储为gm1.m %x = [5999,5903,5848,5700,7884];gm1(x); 测试数据 %二次拟合预测GM(1,1)模型 function gmcal=gm1(x) sizexd2 = size(x,2); %求数组长度 k=0; for y1=x k=k+1; if k>1 x1(k)=x1(k-1)+x(k); %累加生成 z1(k-1)=-0.5*(x1(k)+x1(k-1)); %z1维数减1，用于计算B yn1(k-1)=x(k); else x1(k)=x(k); end end %x1,z1,k,yn1 sizez1=size(z1,2); %size(yn1); z2 = z1'; z3 = ones(1,sizez1)'; YN = yn1'; %转置 %YN B=[z2 z3]; au0=inv(B'*B)*B'*YN; au = au0'; %B,au0,au

ufor = au(2); ua = au(2)./au(1); %afor,ufor,ua %输出预测的 a u 和 u/a的值 constant1 = x(1)-ua; afor1 = -afor; x1t1 = 'x1(t+1)'; estr = 'exp'; tstr = 't'; leftbra = '('; rightbra = ')'; %constant1,afor1,x1t1,estr,tstr,leftbra,rightbra strcat(x1t1,'=',num2str(constant1),estr,leftbra,num2str(afor1),tstr,r ightbra,'+',leftbra,num2str(ua),rightbra) %输出时间响应方程 %****************************************************** %二次拟合 k2 = 0; for y2 = x1 k2 = k2 + 1; if k2 > k else ze1(k2) = exp(-(k2-1)*afor); end end %ze1 sizeze1 = size(ze1,2); z4 = ones(1,sizeze1)'; G=[ze1' z4]; X1 = x1'; au20=inv(G'*G)*G'*X1; au2 = au20'; %z4,X1,G,au20

连续时间混沌系统MATLAB程序和SIMULINK模型

第6章连续时间混沌系统本章讨论连续时间混沌系统的基本特点与分析方法，主要包括混沌数值仿真和硬件实验方法简介、混沌系数平衡点的计算、平衡点的分类与性质、相空间中的轨道、几类典型连续混沌系统的介绍、混沌机理的分析方法、用特征向量空间法寻找异宿轨道、Lorenz系统及混沌机理定性分析、Lorenz映射、Poincare截面、Chua系统及其混沌机理定性分析、时间序列与相空间重构等内容。 6.1 混沌数值仿真和硬件实验方法简介混沌的数值仿真主要包括MA TLAB编程、SIMULINK模块构建、EWB仿真以及其他一些相关的软件仿真或数值计算等方法，从而获取混沌吸引子的相图、时域波形图、李氏指数、分叉图和功率谱等。混沌的硬件实验主要包括模拟/数字电路设计与硬件实验、现场可编程门阵列器件（FPGA）、数字信号处理器（DSP）等硬件实现方法来产生混沌信号。本节仅对各种数值仿真方法作简单介绍。 1）混沌系统的MA TLAB数值仿真该方法主要根据混沌系统的状态方程来编写MA TLAB程序。现举二例来说明这种编程方法。（1）已知Lorenz系统的状态方程为 dx/dt=-a(x-y) dy/dt=bx-xz-y dz/dt=-cz+xy 式中a=10,b=30,c=8/3。 MA TLAB仿真程序如下： >> %************************************************** Function dxdt=lorenz(t,x) %除符号dxdt外，还可用其他编程者习惯的有意义的符号 A=10; B=30; C=8/3; dxdt=zeros(3,1); dxdt(1)=-A*(x(1)-x(2)); dxdt(2)=B*x(1)-x(1).*x(3)-x(2); dxdt(3)=x(1)*x(2)-C*x(3); %************************************************* options=odeset('RelTol',1e-6,'AbsTol',[ 1e-6 1e-6 1e-6]); t0=[0 200]; x0=[0.02,0.01,0.03]; [t,x]=ode45('lorenz',t0,x0,options); %************************************************** n=length(t) n1=round(n/2) %n1=1; %************************************************** figure(1); plot(t(n1:n,1),x(n1:n,1));

灰色预测模型matlab程序精确版

%x=[1019,1088,1324,1408,1601];gm1(x); 测试数据%二次拟合预测GM(1,1)模型 function gmcal=gm1(x) if nargin==0 x=[1019,1088,1324,1408,1601] end format long g sizex=length(x); %求数组长度 k=0; for y1=x k=k+1; if k>1 x1(k)=x1(k-1)+x(k); %累加生成 z1(k-1)=-0.5*(x1(k)+x1(k-1)); %z1维数减1，用于计算B yn1(k-1)=x(k); else x1(k)=x(k); end end %x1,z1,k,yn1 sizez1=length(z1); %size(yn1); z2 = z1'; z3 = ones(1,sizez1)'; YN = yn1'; %转置 %YN B=[z2 z3]; au0=inv(B'*B)*B'*YN; au = au0'; %B,au0,au afor = au(1); ufor = au(2); ua = au(2)./au(1); %afor,ufor,ua %输出预测的 a u 和 u/a的值 constant1 = x(1)-ua; afor1 = -afor; x1t1 = 'x1(t+1)'; estr = 'exp'; tstr = 't'; leftbra = '(';

rightbra = ')'; %constant1,afor1,x1t1,estr,tstr,leftbra,rightbra strcat(x1t1,'=',num2str(constant1),estr,leftbra,num2str(afor1),tstr,rightbra,'+ ',leftbra,num2str(ua),rightbra) %输出时间响应方程 %****************************************************** %二次拟合 k2 = 0; for y2 = x1 k2 = k2 + 1; if k2 > k else ze1(k2) = exp(-(k2-1)*afor); end end %ze1 sizeze1=length(ze1); z4 = ones(1,sizeze1)'; G=[ze1' z4]; X1 = x1'; au20=inv(G'*G)*G'*X1; au2 = au20'; %z4,X1,G,au20 Aval = au2(1); Bval = au2(2); %Aval,Bval %输出预测的 A,B的值 strcat(x1t1,'=',num2str(Aval),estr,leftbra,num2str(afor1),tstr,rightbra,'+',lef tbra,num2str(Bval),rightbra) %输出时间响应方程 nfinal = sizex-1 + 1;（其中+1可改为+5等其他数字，即可预测更多的数字） %决定预测的步骤数5 这个步骤可以通过函数传入 %nfinal = sizexd2 - 1 + 1; %预测的步骤数 1 for k3=1:nfinal x3fcast(k3) = constant1*exp(afor1*k3)+ua; end %x3fcast %一次拟合累加值 for k31=nfinal:-1:0 if k31>1 x31fcast(k31+1) = x3fcast(k31)-x3fcast(k31-1); else if k31>0

Matlab实现混沌系统的控制

基于MATLAB 的各类混沌系统的计算机模拟混沌是非线性系统所独有且广泛存在的一种非周期运动形式, 其覆盖面涉及到自然科学和社会科学的几乎每一个分支。1972年12月29日，美国麻省理工学院教授、混沌学开创人之一E.N.洛伦兹在美国科学发展学会第139次会议上发表了题为《蝴蝶效应》的论文，提出一个貌似荒谬的论断：在巴西一只蝴蝶翅膀的拍打能在美国得克萨斯州产生一个龙卷风，并由此提出了天气的不可准确预报性。为什么会出现这种情况呢？这是混沌在作怪！ “混沌”译自英语中“chaos”一词，原意是混乱、无序，在现代非线性理论中，混沌则是泛指在确定体系中出现的貌似无规则的、类随机的运动。混沌现象是普遍的，就在我们身边，是与我们关系最密切的现象，我们就生活在混沌的海洋中。一支燃着的香烟，在平稳的气流中缓缓升起一缕青烟，突然卷成一团团剧烈搅动的烟雾，向四方飘散；打开水龙头，先是平稳的层流，然后水花四溅，流动变的不规则，这就是湍流；一个风和日丽的夏天，突然风起云涌，来了一场暴风雨。一面旗帜在风中飘扬，一片秋叶从树上落下，它们都在做混沌运动。可见混沌始终围绕在我们的周围，一直与人类为伴。 1.混沌的基本概念 1. 混沌: 目前尚无通用的严格的定义, 一般认为,将不是由随机性外因引起的, 而是由确定性方程(内因)直接得到的具有随机性的运动状态称为混沌。 2. 相空间: 在连续动力系统中, 用一组一阶微分方程描述运动, 以状态变量(或状态向量)为坐标轴的空间构成系统的相空间。系统的一个状态用相空间的一个点表示, 通过该点有唯一的一条积分曲线。 3. 混沌运动: 是确定性系统中局限于有限相空间的高度不稳定的运动。所谓轨道高度不稳定, 是指近邻的轨道随时间的发展会指数地分离。由于这种不稳定性, 系统的长时间行为会显示出某种混乱性。 4. 分形和分维: 分形是 n 维空间一个点集的一种几何性质, 该点集具有无限精细的结构, 在任何尺度下都有自相似部分和整体相似性质, 具有小于所在空间维数 n 的非整数维数。分维就是用非整数维——分数维来定量地描述分形的基本性质。 5. 不动点: 又称平衡点、定态。不动点是系统状态变量所取的一组值, 对于这些值系统不随时间变化。在连续动力学系统中, 相空间中有一个点0x , 若满足当 t →∞时, 轨迹0()x t x →, 则称0x 为不动点。 6. 吸引子: 指相空间的这样的一个点集 s (或一个子空间) , 对s 邻域的几乎任意一点, 当t →∞时所有轨迹线均趋于s, 吸引子是稳定的不动点。 7. 奇异吸引子: 又称混沌吸引子, 指相空间中具有分数维的吸引子的集合。该吸引集由永不重复自身的一系列点组成, 并且无论如何也不表现出任何周期性。混沌轨道就运行在其吸引子集中。 8. 分叉和分叉点: 又称分岔或分支。指在某个或者某组参数发生变化时, 长时间动力学运动的类型也发生变化。这个参数值(或这组参数值)称为分叉点, 在分叉点处参数的微小变化会产生不同性质的动力学特性, 故系统在分叉点处是结构不稳定的。 9. 周期解: 对于系统1()n n x f x += , 当n →∞时,若存在n i n x x ξ+== , 则称该系统有周期i 解ξ 。不动点可以看作是周期为1的解, 因为它满足1n n x x +=。 10. 初值敏感性：对初始条件的敏感依赖是混沌的基本特征，也有人用它来定义混沌：混沌系统是其终极状态极端敏感地依赖于系统的初始状态的系统。敏感依赖性的一个严重后果就在于，使得系统的长期行为变得不可预见。

灰色预测MATLAB程序

灰色预测专设工⑼他QA—叫吋)为原始数列.其1次累 ?加生成数列为恥=妙①曲⑵，…卅何),其中 X° 仇)二工* ° (0.址=1=2= -:n 5-1 卷定义卫的决导数为 d(k) = *町(上)=x 叫咼-x cl)(Jt-l). 令为数列工①的邻值生成数列.即却(去)=^(*) + (1- a)x0)(t-lX 于是定义GM (L 1)的灰微分方程模型为 d(k)-血⑴住)=K 即或严＞(￡) + “尹⑻=人⑴ 在式(1)中』。＞(灼称为灰导数，我称为发展系数, 弧称为白化背景值，b称为灰作用量乜将时刻表殳二2「3「/代入(1)式有 V!1「— a y=代⑶ B =I b*- : X闵0)-Z，：］(K)1

于是G\I <1?1）複至可表示为Y = Bu. 現在问题归结为求sb 在值。用一元线性回归?即最小二秦法求它们的活计值为注二实陌上回归分析中求估计值是用软件计尊的?有标准程序求解，iOmaClab 等。 GM <1? 1>的白化晏対于G\I <1> 1）的灰微分方程（1） >如果将灰导数打（Q 的时刻视为连绫变里"则x°）视为时问（函数卅⑺,于是*〉（Q 対血于导数里级心2 >白化背臬值申的对应于导数卅⑴。于是G\I （1，1）的坝徽分方樂対应于的白微分方程为内?则数堀列X?可以塗互G\I <19 1）且可以进行页色预测。否朋，対数摄做适当的克换处理■如平移叢换: 取C 使得鞍据列严伙）=工⑴伙）+ G 上=1,2,…, 的级比都華住可吝禎盖内。心⑴⑴ + o?i> (r)二dr <2) GM mi ）质色预测的步骤 1 ?教摇的枪绘与处連为了ftilGAl （1,1）建複方法的可行性，亲要为已知期S 做必要的检蛉处理。设原始教据列为了逛=（乂°（1）*6（2）严炉00； >计算数列的级比如果所有的级比都落在可容覆盖区间 ? fc = A- 2,3"?

用Matlab观察分岔与混沌现象

结果分析：在λ取值在[0,0.3]区间内时，y的值保持在0，然后开始上升，在λ取值在0.75附近时，开始分岔为两支。从整体上看，随着λ的值越来越大，所产生的迭代序列越来越复杂，可能会随机地落在区间（-3,3）的任一子区间内。并可能重复，这就是混沌的遍历性。进一步分析：由于λ的取值空间偏小，考虑扩大其取值范围到[0,6]，再进一步观察图像。程序代码如下： clear;clf; hold on axis([0,6,-6,6]); grid for a=0:0.05:6 x=[0.1]; for i=2:150 x(i)=a*sin(pi*x(i-1)); end

灰色预测MATLAB程序

灰色预测M A T L A B程序文档编制序号：[KKIDT-LLE0828-LLETD298-POI08]

灰色预测作用：求累加数列、求a b的值、求预测方程、求残差 clc %清屏，以使结果独立显示 x=[ ]; format long; %设置计算精度 if length(x(:,1))==1 %对输入矩阵进行判断，如不是一维列矩阵，进行转置变换 x=x'; end n=length(x); %取输入数据的样本量 z=0; for i=1:n %计算累加值，并将值赋予矩阵be z=z+x(i,:); be(i,:)=z; end for i=2:n %对原始数列平行移位 y(i-1,:)=x(i,:); end for i=1:n-1 %计算数据矩阵B的第一列数据 c(i,:)=*(be(i,:)+be(i+1,:)); end for j=1:n-1 %计算数据矩阵B的第二列数据 e(j,:)=1; end for i=1:n-1 %构造数据矩阵B B(i,1)=c(i,:); B(i,2)=e(i,:); end alpha=inv(B'*B)*B'*y; %计算参数矩阵即a b的值 for i=1:n+1 %计算数据估计值的累加数列，如改为n+1为n+m可预测后m-1个值 ago(i,:)=(x(1,:)-alpha(2,:)/alpha(1,:))*exp(-alpha(1,:)*(i- 1))+alpha(2,:)/alpha(1,:);%显示输出预测值的累加数列 end var(1,:)=ago(1,:) %显示输出预测值 for i=1:n %如改n为n+m-1，可预测后m-1个值 var(i+1,:)=ago(i+1,:)-ago(i,:); %估计值的累加数列的还原，并计算出下一预测值end for i=1:n error(i,:)=x(i,:)-var(i,:); %计算残差 end c=std(error)/std(x); %调用统计工具箱的标准差函数计算后验差的比值c ago %显示输出预测值的累加数列

灰色预测模型的MATLAB 程序及检验程序

灰色预测模型的Matlab程序及检验程序%灰色预测模型程序 clear syms a b; c=[a b]'; A=[46.232.626.723.020.018.917.516.3];%原始序列B=cumsum(A);%累加 n=length(A); for i=1:(n-1) C(i)=(B(i)+B(i+1))/2; end %计算待定参数 D=A; D(1)=[]; D=D'; E=[-C;ones(1,n-1)]; c=inv(E*E')*E*D; c=c'; a=c(1); b=c(2); %预测往后预测5个数据 F=[];F(1)=A(1); for i=2:(n+5) F(i)=(A(1)-b/a)/exp(a*(i-1))+b/a; end G=[];G(1)=A(1); for i=2:(n+5) G(i)=F(i)-F(i-1); end t1=2002:2009; t2=2002:2014; G plot(t1,A,'o',t2,G) %灰色预测模型检验程序 function[q,c,p]=checkgm(x0,x1) %GM检验函数 %x0原始序列 %x1预测序列 %·返回值

%q–-相对误差 %c--·方差比 %p--小误差概率 e0=x0-x1; q=e0/x0; s1=var(x0); %qpa=mean(e0); s2=var(e0); c=s2/s1; len=length(e0); p=0; for i=1:len if(abs(e0(i))<0.6745*s1) p=p+1; end end p=p/len; end 等级相对误差q方差比C小误差概论P I级<0.01<0.35>0.95 II级<0.05<0.50<0.80 III级<0.10<0.65<0.70 IV级>0.20>0.80<0.60

灰色预测模型matlab程序精确版

灰色预测模型matlab程序灰色模型预测是在数据不呈现一定规律下可以采取的一种建模和预测方法，其预测数据与原始数据存在一定的规律相似性 %下面程序是灰色模型GM(1,1)程序二次拟合和等维新陈代谢改进预测程序,mat lab6.5 ,使用本程序请注明，程序存储为gm1.m %x = [5999,5903,5848,5700,7884];gm1(x); 测试数据 %二次拟合预测GM(1,1)模型 function gmcal=gm1(x) sizexd2 = size(x,2); %求数组长度 k=0; for y1=x k=k+1; if k>1 x1(k)=x1(k-1)+x(k); %累加生成 z1(k-1)=-0.5*(x1(k)+x1(k-1)); %z1维数减1，用于计算B yn1(k-1)=x(k); else x1(k)=x(k); end end %x1,z1,k,yn1 sizez1=size(z1,2); %size(yn1); z2 = z1'; z3 = ones(1,sizez1)'; YN = yn1'; %转置 %YN B=[z2 z3]; au0=inv(B'*B)*B'*YN;

%B,au0,au afor = au(1); ufor = au(2); ua = au(2)./au(1); %afor,ufor,ua %输出预测的 a u 和 u/a的值 constant1 = x(1)-ua; afor1 = -afor; x1t1 = 'x1(t+1)'; estr = 'exp'; tstr = 't'; leftbra = '('; rightbra = ')'; %constant1,afor1,x1t1,estr,tstr,leftbra,rightbra strcat(x1t1,'=',num2str(constant1),estr,leftbra,num2str(afor1),tstr,r ightbra,'+',leftbra,num2str(ua),rightbra) %输出时间响应方程 %****************************************************** %二次拟合 k2 = 0; for y2 = x1 k2 = k2 + 1; if k2 > k else ze1(k2) = exp(-(k2-1)*afor); end end %ze1 sizeze1 = size(ze1,2); z4 = ones(1,sizeze1)'; G=[ze1' z4]; X1 = x1'; au20=inv(G'*G)*G'*X1;

数字分析 matlab程序

subplot(2,2,1) x=0.578;c=0.5; hold on; for k=1:60 x y=c*x*(1-x); y plot(x,y,'bd') title('研究一般迭代公式的复杂行为混沌现象') x=y; k=k+1 end subplot(2,2,2) x=0.578;c=1.5; hold on; for k=1:60 x y=c*x*(1-x); y plot(x,y,'g*') title('研究一般迭代公式的复杂行为混沌现象') x=y; k=k+1 end subplot(2,2,3) x=0.578;c=2.5; hold on; for k=1:60 x y=c*x*(1-x); y plot(x,y,'k+') title('研究一般迭代公式的复杂行为混沌现象') x=y; k=k+1 end subplot(2,2,4) x=0.578;c=4; hold on; for k=1:60 x y=c*x*(1-x); y plot(x,y,'ro')

title('研究一般迭代公式的复杂行为混沌现象') x=y; k=k+1 end x = 0.5780 y = 0.1220 k = 2 x = 0.1220 y = 0.0535 k = 3 x = 0.0535 y = 0.0253

k = 4 x = 0.0253 y = 0.0123 k = 5 x = 0.0123 y = 0.0061 k = 6 x = 0.0061 y = 0.0030

k = 7 x = 0.0030 y = 0.0015 k = 8 x = 0.0015 y = 7.5413e-004 k = 9 x = 7.5413e-004 y =

灰色理论预测模型及GM(1,1)matlab程序

灰色理论预测模型及GM(1,1)matlab程序灰色预测方法简介灰色预测是一种对含有不确定因素的系统进行预测的方法。灰色预测通过鉴别系统因素之间发展趋势的相异程度，即进行关联分析，并对原始数据进行生成处理来寻找系统变动的规律，生成有较强规律性的数据序列，然后建立相应的微分方程模型，从而预测事物未来发展趋势的状况。其用等时距观测到的反应预测对象特征的一系列数量值构造灰色预测模型，预测未来某一时刻的特征量，或达到某一特征量的时间。通过对原始数据的整理寻找数的规律，分为三类： a、累加生成：通过数列间各时刻数据的依个累加得到新的数据与数列。累加前数列为原始数列，累加后为生成数列。 b、累减生成：前后两个数据之差，累加生成的逆运算。累减生成可将累加生成还原成非生成数列。 c、映射生成：累加、累减以外的生成方式。建模步骤 a、建模机理 b、把原始数据加工成生成数； c、对残差（模型计算值与实际值之差）修订后，建立差分微分方程模型； d、基于关联度收敛的分析； e、gm模型所得数据须经过逆生成还原后才能用。 f、采用“五步建模（系统定性分析、因素分析、初步量化、动态量化、优化）”法，建立一种差分微分方程模型gm(1,1）预测模型。 GM(1,1)程序： % 本程序主要用来计算根据灰色理论建立的模型的预测值。 % 应用的数学模型是GM(1,1)。 % 原始数据的处理方法是一次累加法。 clear;clc; % load ('data.txt');

% y=data'; y=[3 4 5 4 7 7]; n=length(y); yy=ones(n,1); yy(1)=y(1); for i=2:n yy(i)=yy(i-1)+y(i); end B=ones(n-1,2); for i=1:(n-1) B(i,1)=-(yy(i)+yy(i+1))/2; B(i,2)=1; end BT=B'; for j=1:n-1 YN(j)=y(j+1); end YN=YN'; A=inv(BT*B)*BT*YN; a=A(1); u=A(2); t=u/a; t_test=input('请输入需要预测个数：'); i=1:t_test+n; yys(i+1)=(y(1)-t).*exp(-a.*i)+t; yys(1)=y(1); for j=n+t_test:-1:2 ys(j)=yys(j)-yys(j-1); end x=1:n; xs=2:n+t_test; yn=ys(2:n+t_test); plot(x,y,'^r',xs,yn,'*-b'); det=0; for i=2:n det=det+abs(yn(i)-y(i)); end det=det/(n-1); disp(['百分绝对误差为：',num2str(det),'%']); disp(['预测值为：',num2str(ys(n+1:n+t_test))]);

灰色系统预测GM(1-1)模型及其Matlab实现

灰色系统预测GM(1,1)模型及其Matlab 实现预备知识 (1)灰色系统白色系统是指系统内部特征是完全已知的；黑色系统是指系统内部信息完全未知的；而灰色系统是介于白色系统和黑色系统之间的一种系统，灰色系统其内部一部分信息已知，另一部分信息未知或不确定。 (2)灰色预测灰色预测，是指对系统行为特征值的发展变化进行的预测，对既含有已知信息又含有不确定信息的系统进行的预测，也就是对在一定范围内变化的、与时间序列有关的灰过程进行预测。尽管灰过程中所显示的现象是随机的、杂乱无章的，但毕竟是有序的、有界的，因此得到的数据集合具备潜在的规律。灰色预测是利用这种规律建立灰色模型对灰色系统进行预测。目前使用最广泛的灰色预测模型就是关于数列预测的一个变量、一阶微分的GM(1,1)模型。它是基于随机的原始时间序列，经按时间累加后所形成的新的时间序列呈现的规律可用一阶线性微分方程的解来逼近。经证明，经一阶线性微分方程的解逼近所揭示的原始时间序列呈指数变化规律。因此，当原始时间序列隐含着指数变化规律时，灰色模型GM(1,1)的预测是非常成功的。 1 灰色系统的模型GM(1,1) 1.1 GM(1,1)的一般形式设有变量X (0)＝{X (0)(i)，i=1,2，...，n}为某一预测对象的非负单调原始数据列，为建立灰色预测模型：首先对X (0)进行一次累加(1—AGO, Acumulated Generating Operator)生成一次累加序列： X (1)＝{X (1)(k )，k ＝1，2，…，n} 其中 X (1)(k )＝ ∑ =k i 1 X (0)(i) ＝X (1)(k －1)+ X (0)(k ) (1) 对X (1)可建立下述白化形式的微分方程： dt dX )1(十) 1(aX ＝u (2) 即GM(1,1)模型。上述白化微分方程的解为(离散响应)： ∧ X (1)(k +1)＝(X (0)(1)－ a u )ak e -＋a u (3) 或 ∧ X (1)(k )＝(X (0)(1)－ a u ))1(--k a e ＋a u (4)