M文件和Simulink求解连续微分系统实例分析

matlab simulink求解微分方程组例题以下是一个使用MATLAB Simulink求解微分方程组的例子：假设有一个简单的二阶线性微分方程组：dx1/dt = a11*x1 + a12*x2dx2/dt = a21*x1 + a22*x2其中，a11、a12、a21和a22是给定的常数。

1. 创建一个新的Simulink模型。

2. 在模型中添加两个积分器（Integrator）模块，分别表示x1和x2的积分变量。

3. 添加两个乘法器（Gain）模块，分别表示a11和a12，将常数a11和a12连接到相应的乘法器输入。

4. 将积分器和乘法器的输出连接到一个加法器（Sum）模块，表示a11*x1 + a12*x2。

5. 重复步骤3和4，添加两个乘法器和一个加法器，表示a21*x1 + a22*x2。

6. 将两个加法器的输出连接到相应的积分器的输入，表示dx1/dt和dx2/dt。

7. 添加一个作用于积分器输出的数据记录器（To Workspace）模块，用于记录x1和x2的值。

通过上述步骤，我们建立了一个包含两个积分器和四个乘法器的模型，表示二阶线性微分方程组。

我们可以通过设置初始条件和常数a11、a12、a21和a22的值，使用Simulink模型求解微分方程组并观察结果。

例如，假设初始条件为x1(0) = 1，x2(0) = 0，a11 = 2，a12 = -1，a21 = 1，a22 = -2。

设置初始条件后，运行Simulink模型，模型将在仿真时间内求解微分方程组，并将x1和x2的值保存在工作空间中。

通过查看工作空间中的数据记录，我们可以获得仿真时间内x1和x2的值，并观察它们的变化情况。

Simulink是MATLAB的一个模块，主要用于建立、模拟和分析动态系统。

对于微分方程组，Simulink提供了一些工具和模块，如`dsolve`和`ode45`等，可以帮助我们求解微分方程组。

以下是一个使用Simulink求解微分方程组的步骤：
1. **建立微分方程组**：首先，我们需要定义一个微分方程组。

例如，考虑以下微分方程组：
dx/dt = y
dy/dt = -x
2. **使用dsolve求解**：在MATLAB命令窗口中，我们可以使用`dsolve`函数来求解这个微分方程组。

例如：
`syms x y t
dsolve([x' == y, y' == -x], [x, y])`
这将返回微分方程组的解析解。

3. **在Simulink中建模**：然后，我们可以在Simulink中建立模型。

首先，我们需要创建一个新模型，然后在模型中添加`Scope`模块来显示结果。

4. **添加微分方程求解器**：在Simulink的库浏览器中，找到并添加`Math Operations` -> `Differential Equations` -> `ode45`模块。

这个模块可以用来求解微分方程。

5. **设置参数**：双击`ode45`模块，设置参数。

在`Function`选项中，选择我们之前定义的微分方程组。

在`Initial conditions`选项中，设置初始条件。

6. **运行模型**：最后，运行模型，观察结果。

注意：以上步骤是一个简化的示例，实际操作可能需要根据具体问题进行调整。

MATLAB Simulink是MATLAB的一个强大的图形用户界面，可用于设计和模拟动态系统。

在Simulink中，可以使用不同的工具和库来求解微分方程，包括但不限于：1. **Simulink Coder**：这是一个用于生成嵌入式代码的强大工具，可以用于解决微分方程。

2. **Simulink Solver**：这是Simulink自带的一个求解器，可以用于解决常微分方程。

3. **Simulink Simulation**：这是Simulink的基本求解器，可以用于解决各种动态系统问题，包括微分方程。

以下是一个简单的例子，说明如何在Simulink中解决微分方程：1. **打开Simulink**：在MATLAB中，点击“新建”->“Simulink 模型”。

2. **创建模型**：在Simulink编辑器中，你可以添加各种模块来构建你的模型。

在这个例子中，我们首先需要一个模块来表示我们的微分方程。

我们可以使用`Simulink.Continuous`模块来表示一个连续系统，然后在该模块的参数设置中输入我们的微分方程。

例如，如果你的微分方程是 `dy/dx = f(x, y)`，你可以在参数设置中输入这个方程。

3. **设置初始条件和边界条件**：你可能还需要设置初始条件和边界条件。

你可以通过添加`InitialCondition`和`BoundaryCondition`模块来实现这一点。

4. **设置求解器**：在Simulink编辑器的顶部菜单中，选择“模型配置参数”。

在弹出的窗口中，选择“求解器”选项卡，然后选择你需要的求解器。

例如，你可以选择“auto”让Simulink自动选择合适的求解器，或者你可以选择“ode45”来使用MATLAB内置的求解器。

5. **运行模拟**：一旦你设置了所有必要的参数和模块，你就可以运行模拟了。

在Simulink编辑器的顶部菜单中，选择“模拟”->“开始”。

实验二 数字PID 控制计算机控制是一种采样控制,它只能根据采样时刻的偏差值计算控制量;因此连续PID 控制算法不能直接使用,需要采用离散化方法;在计算机PID 控制中,使用的是数字PID 控制器;一、位置式PID 控制算法按模拟PID 控制算法,以一系列的采样时刻点kT 代表连续时间t,以矩形法数值积分近似代替积分,以一阶后向差分近似代替微分,可得离散PID 位置式表达式：式中,D p d I p i T k k T k k ==,,e 为误差信号即PID 控制器的输入,u 为控制信号即控制器的输出;在仿真过程中,可根据实际情况,对控制器的输出进行限幅;二、连续系统的数字PID 控制仿真连续系统的数字PID 控制可实现D/A 及A/D 的功能,符合数字实时控制的真实情况,计算机及DSP 的实时PID 控制都属于这种情况;1．Ex3 设被控对象为一个电机模型传递函数Bs Js s G +=21)(,式中J=,B=;输入信号为)2sin(5.0t π,采用PD 控制,其中5.0,20==d p k k ;采用ODE45方法求解连续被控对象方程;因为Bs Js s U s Y s G +==21)()()(,所以u dt dy B dty d J =+22,另y y y y ==2,1,则⎪⎩⎪⎨⎧+-==/J)*u ((B/J)y y y y 12221 ,因此连续对象微分方程函数如下 function dy = ex3ft,y,flag,parau=para;J=;B=;dy=zeros2,1;dy1 = y2;dy2 = -B/Jy2 + 1/Ju;控制主程序clear all;close all;ts=; %采样周期xk=zeros2,1;%被控对象经A/D转换器的输出信号y的初值e_1=0;%误差ek-1初值u_1=0;%控制信号uk-1初值for k=1:1:2000 %k为采样步数timek = kts; %time中存放着各采样时刻rink=sin12pikts; %计算输入信号的采样值para=u_1; % D/AtSpan=0 ts;tt,xx=ode45'ex3f',tSpan,xk,,para; %ode45解系统微分方程%xx有两列,第一列为tt时刻对应的y,第二列为tt时刻对应的y导数xk = xxend,:; % A/D,提取xx中最后一行的值,即当前y和y导数youtk=xk1; %xk1即为当前系统输出采样值ykek=rink-youtk;%计算当前误差dek=ek-e_1/ts; %计算uk中微分项输出uk=ek+dek;%计算当前uk的输出%控制信号限幅if uk>uk=;endif uk<uk=;end%更新uk-1和ek-1u_1=uk;e_1=ek;endfigure1;plottime,rin,'r',time,yout,'b';%输入输出信号图xlabel'times',ylabel'rin,yout';figure2;plottime,rin-yout,'r';xlabel'times',ylabel'error';%误差图程序运行结果显示表1所示;表1 程序运行结果分析:输出跟随输入,PD 控制中,微分控制可以改善动态特性,调节时间缩短,允许加大比例控制,使稳态误差减小,提高了控制精度.2．Ex4 被控对象是一个三阶传递函数ss s 1047035.8752350023++,采用Simulink 与m 文件相结合的形式,利用ODE45方法求解连续对象方程,主程序由Simulink 模块实现,控制器由m 文件实现;输入信号为一个采样周期1ms 的正弦信号;采用PID 方法设计控制器,其中05.0,2,5.1===d i p k k k ;误差初始化由时钟功能实现,从而在m 文件中实现了误差的积分和微分;控制主程序：控制子程序：function u=ex4fu1,u2%u1为Clock,u2为图2-1中Sum 模块输出的误差信号e 的采样值 persistent errori error_1if u1==0 %当Clock=0时,即初始时,ek=ek-1=0errori=0error_1=0endts=;kp=;ki=;kd=;error=u2;errord=error-error_1/ts;%一阶后向差分误差信号表示的误差微分errori=errori+errorts;%累积矩形求和计算的误差的积分u=kperror+kderrord+kierrori;%由PID 算式得出的当前控制信号ukerror_1=error;%误差信号更新图2-1 Simulink 仿真程序其程序运行结果如表2所示;Matlab 输出结果errori =error_1 =表2 例4程序运行结果三、离散系统的数字PID 控制仿真1．Ex5 设被控对象为ss s s G 1047035.87523500)(23++=,采样时间为1ms,对其进行离散化;针对离散系统的阶跃信号、正弦信号和方波信号的位置响应,设计离散PID 控制器;其中S 为信号选择变量,S=1时是阶跃跟踪,S=2时为方波跟踪,S=3时为正弦跟踪;求出Gs 对应的离散形式)()()(z U z Y z G =,其中Yz 和Uz 是关于z 的多项式,则可以得到其对应的差分表达式仿真程序：%PID Controllerclear all;close all;ts=;%采样周期sys=tf,1,,,0;%被控对象连续传递函数dsys=c2dsys,ts,'z';%转换成离散z 传递函数的形式num,den=tfdatadsys,'v';%提取z 传递函数中的分子和分母多项式系数u_1=;u_2=;u_3=;%uk-1、uk-2、uk-3的初值y_1=;y_2=;y_3=; %yk-1、yk-2、yk-3的初值x=0,0,0'; %比例、微分、积分项的初值error_1=0;%ek-1的初值disp'S=1--step,S=2--sin,S=3--square'% S=1阶跃,S=2方波,S=3正弦S=input'Number of input signal S:'%接收输入信号代号for k=1:1:1500timek=kts;%各采样时刻if S==1 %阶跃输入时kp=;ki=;kd=; %各项PID系数rink=1; %阶跃信号输入elseif S==2kp=;ki=;kd=; %各项PID系数rink=signsin22pikts; %方波信号输入elseif S==3kp=;ki=;kd=; %各项PID系数rink=sin22pikts; %正弦信号输入enduk=kpx1+kdx2+kix3; %PID控制信号输出uk%控制信号输出限幅if uk>=10uk=10;endif uk<=-10uk=-10;end%根据差分方程计算系统当前输出ykyoutk=-den2y_1-den3y_2-den4y_3+num2u_1+num3u_2+num4u_3; errork=rink-youtk;%当前误差%更新uk-1、uk-2、uk-3、yk-1、yk-2、yk-3u_3=u_2;u_2=u_1;u_1=uk;y_3=y_2;y_2=y_1;y_1=youtk;x1=errork; %比例输出x2=errork-error_1/ts; %微分输出x3=x3+errorkts; %积分输出error_1=errork; %更新ek-1endfigure1; %作图plottime,rin,'r',time,yout,'b';xlabel'times',ylabel'rin,yout';其程序运行结果如表3所示;2．Ex6针对于Ex5被控对象所对应的离散系统,设计针对三角波、锯齿波和随机信号的位置式响应;仿真程序：;程序中当S=1时为三角波,S=2时为锯齿波,S=3时为随机信号;如果D=1,则通过pause命令实现动态演示仿真;%PID Controllerclear all;close all;ts=;sys=tf,1,,,0;dsys=c2dsys,ts,'z';num,den=tfdatadsys,'v';u_1=;u_2=;u_3=;r_1=rand;y_1=0;y_2=0;y_3=0;x=0,0,0';error_1=0;disp'S=1--Triangle,S=2--Sawtooth,S=3--Random'% S=1三角,S=2锯齿,S=3随机S=input'Number of input signal S:'%接收输入信号代号disp'D=1--Dynamic display,D~=1--Direct display'%D=1动画显示,D~=1直接显示D=input'D='for k=1:1:3000timek=kts;kp=;ki=;kd=;if S==1 %Triangle Signalif modtimek,2<1rink=modtimek,1;elserink=1-modtimek,1;endrink=rink;endif S==2 %Sawtooth Signalrink=modtimek,;endif S==3 %Random Signalrink=rand;vrk=rink-r_1/ts; %Max speed iswhile absvrk>=rink=rand;vrk=absrink-r_1/ts;endenduk=kpx1+kdx2+kix3; %PID Controller%Restricting the output of controllerif uk>=10uk=10;endif uk<=-10uk=-10;end%Linear modelyoutk=-den2y_1-den3y_2-den4y_3+num2u_1+num3u_2+num4u_3; errork=rink-youtk;r_1=rink;u_3=u_2;u_2=u_1;u_1=uk;y_3=y_2;y_2=y_1;y_1=youtk;x1=errork; %Calculating Px2=errork-error_1/ts; %Calculating Dx3=x3+errorkts; %Calculating Ixik=x3;error_1=errork;if D==1 %Dynamic Simulation Displayplottime,rin,'b',time,yout,'r';pause;endendplottime,rin,'r',time,yout,'b';xlabel'times';ylabel'rin,yout';Matlab运行结果为：S=1--Triangle,S=2--Sawtooth,S=3--RandomNumber of input signal S:12、3S =1D=1--Dynamic display,D=0--Direct display %D=1动画显示,D~=1直接显示D=0D =0 %D=0直接显示,如果D=1,则通过pause命令实现动态演示仿真;其程序运行结果如表4所示;表4 程序运行结果分析：s=1时为三角波的位置式响应,s=2时为锯齿波的位置式响应；s=3时为随机信号的位置式响应;3．Ex7 采用Simulink实现PID控制器的设计,如图2-2所示,其中离散PID控制的子系统如图2-3所示,其封装界面如图2-4所示;仿真程序：图2-2 离散PID控制的Simulink主程序图2-3 离散PID控制的Simulink控制器程序图2-4 离散PID控制的封装界面其程序运行结果如表5所示;分析：位置式PID控制算法的缺点是,由于采用全量输出,所以每次输出均与过去的状态有关,计算时要对ek量进行累加,计算机输出控制量uk对应的是执行机构的实际位置偏差,如果位置传感器出现故障,uk 可能会出现大幅度变化;uk 大幅度变化会引起执行机构未知的大幅度变化,这种情况在生产中是不允许的,在某些重要场合还可能造成重大事故;为了避免这种情况的发生,可采用增量式PID 控制算法;四、增量式PID 控制算法及仿真当执行机构需要的是控制量的增量例如驱动步进电机时,应采用增量式PID 控制,根据递推原理可得增量式PID 控制算法为Ex8 设被控对象ss s G 50400)(2+=,采用增量式控制算法,PID 控制参数10,1.0,8===d i p k k k ; 仿真程序：%Increment PID Controllerclear all;close all;ts=;sys=tf400,1,50,0;dsys=c2dsys,ts,'z';num,den=tfdatadsys,'v';u_1=;u_2=;u_3=;y_1=0;y_2=0;y_3=0;x=0,0,0';error_1=0;error_2=0;for k=1:1:1000timek=kts;rink=;kp=8;ki=;kd=10;duk=kpx1+kdx2+kix3;uk=u_1+duk;if uk>=10 uk=10; end if uk<=-10 uk=-10; endyoutk=-den2y_1-den3y_2+num2u_1+num3u_2; error=rink-youtk;u_3=u_2;u_2=u_1;u_1=uk; y_3=y_2;y_2=y_1;y_1=youtk;x1=error-error_1; %Calculating P x2=error-2error_1+error_2; %Calculating D x3=error; %Calculating I error_2=error_1; error_1=error; endplottime,rin,'b',time,yout,'r'; xlabel'times';ylabel'rin,yout'; 其程序运行结果如图5所示;分析：增量式PID 控制算法计算的误差小;由于他只输出控制量,所以误动作时影响小,系统的超调及振动也少;由于控制算法中不需要累加,控制增量Δuk 仅与最近k 次的采样有关,所以误动作影响小,而且较容易通过加权处理获得比较好的控制效果;仿真程序运行结果：PID 控制参数10,1.0,8===d i p k k k ;图5 增量式PID 控制算法仿真结果实验三 PID 控制的改进算法在计算机控制系统中,PID 控制是通过计算机程序实现的,因此灵活性很大;一些原来在模拟PID 控制器中无法实现的问题,在引入计算机以后,就可以得到解决,于是产生了一系列的改进算法,形成非标准的控制算法,以改善系统品质,满足不同控制系统的需要;一、积分分离PID 控制算法在普通PID 控制中,积分的目的是为了消除金叉,提高精度,但在过程的启动、结束或大幅度增减设定是,短时间内系统输出有很大偏差,会造成PID 运算的积分积累,致使控制量超过执行机构可能允许的最大动作范围对应的极限控制量,引起系统较大的超调,甚至引起系统较大的振荡,这在生产中是绝对不允许的;积分分离控制基本思路是,当被控量与设定值偏差较大时,取消积分作用,以免由于积分作用使系统稳定性降低,超调量增大；当被控量接近给定值时,引入积分控制,以便消除静差,提高控制精度;其具体实现步骤是：1）根据实际情况,人为设定阈值ε>0；2）当ε>)(k e 时,采用PD 控制,可避免产生过大的超调,又使系统有较快的响应； 3）当ε≤)(k e 时,采用PID 控制,以保证系统的控制精度; 积分分离算法可表示为：式中,T 为采样时间,β为积分项的开关系数,⎩⎨⎧>≤=ξξβ|)(|0|)(|1k e k eEx9 设备控对象为一个延迟对象160)(80+=-s e s G s,采样周期为20s,延迟时间为4个采样周期,即80s;输入信号rk=40,控制器输出限制在-110,110;3,005.0,8.0===d i p k k k 被控对象离散化为)5()2()1()2()(-+--=k u num k y den k y仿真方法一：仿真程序：;当M=1时采用分段积分分离法,M=2时采用普通PID 控制; %Integration Separation PID Controller clear all; close all; ts=20; %Delay plantsys=tf1,60,1,'inputdelay',80; dsys=c2dsys,ts,'zoh'; num,den=tfdatadsys,'v';u_1=0;u_2=0;u_3=0;u_4=0;u_5=0; y_1=0;y_2=0;y_3=0;error_1=0;error_2=0;ei=0;% M=1分段积分分离,M=2普通PIDdisp'M=1--Using integration separation,M=2--Not using integration separation' M=input'whether or not use integration separation method:'for k=1:1:200timek=kts;%输出信号youtk=-den2y_1+num2u_5;rink=40;errork=rink-youtk;ei=ei+errorkts;%积分项输出if M==1 %使用分段积分分离if abserrork>=30&abserrork<=40beta=;elseif abserrork>=20&abserrork<=30beta=;elseif abserrork>=10&abserrork<=20beta=;elsebeta=;endelseif M==2beta=;endkp=;ki=;kd=;uk=kperrork+kderrork-error_1/ts+betakiei;if uk>=110 % 控制信号限幅uk=110; end if uk<=-110 uk=-110; endu_5=u_4;u_4=u_3;u_3=u_2;u_2=u_1;u_1=uk; y_3=y_2;y_2=y_1;y_1=youtk; error_2=error_1; error_1=errork; end figure1;plottime,rin,'b',time,yout,'r'; xlabel'times';ylabel'rin,yout'; figure2; plottime,u,'r'; xlabel'times';ylabel'u'; 其程序运行结果表1所示;表1 例9仿真方法一结果分析：积分饱和的防止方法有积分分离法和预限削弱法;积分作用使系统稳定性降低,超调量增大;比较仿真结果,当被控量与设定值偏差较大时,删除积分作用,以使∑=kj j e 0)(不至过大;只有当)(k e 较小时方引入积分作用,以消除静差,提高控制精度,控制量不宜进入饱和区;仿真方法二：采用Simulink 仿真初始化程序clear all;close all;ts=20;sys=tf1,60,1,'inputdelay',80;dsys=c2dsys,ts,'zoh';num,den=tfdatadsys,'v';kp=;ki=;kd=;Simulink主程序,如图3-1所示;图3-1 Simulink主程序其运行结果如表2所示;表2 Simulink仿真结果分析：由图可知,积分时间常数能消除系统的稳态误差,提高系统控制精度,只有当积分时间常数合适时,过度过程的特性才比较理想;积分时间常数过小,系统震荡次数多,积分时间常数过大,对系统性能影响减少;二、抗积分饱和PID控制算法所谓积分饱和是指若系统存在一个方向的偏差,PID控制器的输出由于积分作用的不断累加而加大,从而导致执行机构达到极限位置Xmax,若控制器输出uk继续增大,阀门开度不可能在增大,此时就称计算机输出控制超出正常运行范围而进入了饱和区;一旦系统出现反向偏差,uk逐渐从饱和区推出;进入饱和区越深,则退出饱和区所需时间越长;在这段时间内,执行机构仍停留在极限位置而不能随偏差反向立即作出相应的改变,这时系统就像失去控制一样,造成控制性能恶化;这种现象称为积分饱和现象或积分失控现象;抗积分饱和的思路是,在计算uk时,首先判断上一时刻的控制量uk-1是否已超出限制范围;若uk-1>u max ,则只累加负偏差；若uk-1<u min ,则只累加正偏差;这种算法可以避免控制量长时间停留在饱和区;Ex10 设被控对象为ss s s G 1047035.87523500)(23++=,采样周期1ms;输入rk=30,0,9,85.0===d i p k k k仿真程序：;M=1时采用抗积分饱和算法,M=2时采用普通PID 算法; %PID Controler with intergration sturation clear all; close all; ts=; sys=tf,1,,,0; dsys=c2dsys,ts,'z'; num,den=tfdatadsys,'v'; u_1=;u_2=;u_3=; y_1=0;y_2=0;y_3=0; x=0,0,0'; error_1=0;um=6;%控制信号限幅值 kp=;ki=;kd=;rin=30; %Step Signal % M=1抗积分饱和,M=2普通PIDdisp'M=1--Using intergration sturation,M=2--Not using iintergration sturation' M=input'whether or not use integration separation method:' for k=1:1:800 timek=kts;uk=kpx1+kdx2+kix3; % PID Controller if uk>=um uk=um; end if uk<=-umuk=-um;end%Linear modelyoutk=-den2y_1-den3y_2-den4y_3+num2u_1+num3u_2+num4u_3; errork=rin-youtk;if M==1 %Using intergration sturationif uk>=umif errork>0alpha=0;elsealpha=1;endelseif uk<=-umif errork>0alpha=1;elsealpha=0;endelsealpha=1;endelseif M==2 %Not using intergration sturationalpha=1;end%Return of PID parametersu_3=u_2;u_2=u_1;u_1=uk;y_3=y_2;y_2=y_1;y_1=youtk;error_1=errork;x1=errork; % 计算比例项x2=errork-error_1/ts; % 计算微分项x3=x3+alphaerrorkts; % 计算积分项 xik=x3; end figure1; subplot311;plottime,rin,'b',time,yout,'r';xlabel'times';ylabel'Position tracking'; subplot312; plottime,u,'r';xlabel'times';ylabel'Controller output'; subplot313; plottime,xi,'r';xlabel'times';ylabel'Integration'; 其运行结果如表3所示;表3 例10仿真结果分析:比较仿真结果知,采用普通的算法时,积分项的存在,有时可能会引起积分饱和,增加系统的调整时间和超调量,而采用了抗积分饱和的方法,可以消除静态误差,使控制量不易进入饱和区,即使进入了,也能较快,系统的输出特性得到了一定改善;三、不完全微分PID 控制算法在PID 控制中,微分信号的引入可改善系统的动态特性,但也易引入高频干扰,在误差扰动突变时尤其显出微分项的不足;若在控制算法中加入低通滤波器,则可使系统性能得到改善;具体做法就是在PID 算法中加入一个一阶惯性环节低通滤波器sT f +11,T f 为滤波器系数;可得此时的微分项输出为()())1()1()()1()1()()1()(-+---=--++-+=k u k e k e K k e k e T T T k k u T T T k u D D f s DpD fs f D αα,其中)1(-+=k u T T T D fs f α,sDpD T T k K =,T s 为采样时间,T D 为微分时间常数; Ex11 被控对象为时滞系统传递函数160)(80+=-s e s G s,在对象的输出端加幅值为的随机信号;采样周期为20ms;采用不完全微分算法,140,0055.0,3.0===D i p T k k ;所加的低通滤波器为11801)(+=s s Q仿真程序：;M=1时采用不完全微分,M=2时采用普通PID %PID Controler with Partial differential clear all; close all; ts=20;sys=tf1,60,1,'inputdelay',80; dsys=c2dsys,ts,'zoh'; num,den=tfdatadsys,'v';u_1=0;u_2=0;u_3=0;u_4=0;u_5=0; %控制信号初值 ud_1=0; %uDk-1初值y_1=0;y_2=0;y_3=0; %输出信号初值 error_1=0; ei=0; for k=1:1:100 timek=kts; rink=;youtk=-den2y_1+num2u_5; %输出信号差分方程 Dk=rands1;%干扰信号youtk=youtk+Dk; %加入干扰后的输出信号 errork=rink-youtk;ei=ei+errorkts; %矩形面积求和计算的积分项输出 kp=; ki=;TD=140;kd=kpTD/ts;Tf=180;%Q的滤波器系数Q=tf1,Tf,1; %低通滤波器%M=1选择不完全微分,M=2选择普通PIDdisp'M=1—Using Partial differential PID,M=2-- Using PID Controler without Partial differential'M=input'whether or not use Partial differential PID:'if M==1 %M=1时用不完全微分alfa=Tf/ts+Tf;udk=kd1-alfaerrork-error_1+alfaud_1;uk=kperrork+udk+kiei;ud_1=udk;elseif M==2 %M=2时用普通PIDuk=kperrork+kderrork-error_1+kiei;end%输出限幅if uk>=10uk=10;endif uk<=-10uk=-10;end%更新采样值u_5=u_4;u_4=u_3;u_3=u_2;u_2=u_1;u_1=uk;y_3=y_2;y_2=y_1;y_1=youtk;error_1=errork;endfigure1;plottime,rin,'b',time,yout,'r';xlabel'times';ylabel'rin,yout'; figure2; plottime,u,'r'; xlabel'times';ylabel'u'; figure3;plottime,rin-yout,'r'; xlabel'times';ylabel'error'; figure4; bodeQ,'r'; dcgainQ;其运行结果如表4所示;分析：比较m=1与m=2的图可得知：在标准PID 算式中,当有阶跃信号输入时,微分想输出急剧增加,容易引起调节过程的震荡,导致品质因数下降,为了克服这点才引入不完全微分的PID 算法;其微分作用逐渐下降,微分输出信号按指数规律逐渐衰减到零,因而系统变化比较缓慢,不容易引起振荡;微分控制可以改善动态特性,如超调量减少,调节时间缩短,使稳态误差减少,提高控制精度;表4 例11运行结果四、微分线性PID 控制算法微分线性的PID 控制结构如图3-2所示,其特点是只对输出量yk 进行微分,而对给定值rk 不进行微分;这样,在改变给定值时,输出不会改变,而被控量的变化通常是比较缓和的,它适用于给定值rk 频繁升降的场合,可以避免给定值升降时引起的系统振荡,从而改善系统的动态特性;图3-2 微分先行PID 控制结构图 令微分部分的传递函数为111)()(<++=γγs T s T s y s u D D D ,式中11+s T D γ相当于低通滤波器;则有y dtdyT u dt du T D D D D+=+γ 由差分得：)()1()()()1()(k y Tk y k y T k u T k u k u T D D D D D+--=+--γ整理得微分部分的输出：)1()()1()(321--+-=k y c k y c k u c k u D D 其中TT T c T T T T c TT T c D DD D D D +=++=+=γγγγ321,,比例积分部分的传递函数为：⎪⎪⎭⎫⎝⎛+=s T k s E s u I p PI 11)()(,其中T I 为积分时间常数; 离散控制算式为)()()(k u k u k u D PI +=;Ex12 设被控对象为一个延迟对象160)(80+=-s e s G s,采样周期为20s;输入信号为带有高频干扰的方波信号：))03.0sin(05.0)0005.0sgn(sin()(t t t r ππ+=;普通PID 控制中 14,0021.0,36.0===d i p k k k ;微分先行PID 中γ=;仿真程序：;M=1时使用微分先行PID,M=2使用普通PID %PID Controler with differential in advance clear all; close all; ts=20;sys=tf1,60,1,'inputdelay',80; dsys=c2dsys,ts,'zoh'; num,den=tfdatadsys,'v';u_1=0;u_2=0;u_3=0;u_4=0;u_5=0; ud_1=0;y_1=0;y_2=0;y_3=0; error_1=0;error_2=0; ei=0;%M=1使用微分先行PID,M=2使用普通PIDdisp'M=1aUsing PID Controler with differential in advance ,M=2-- Using common PID Controler'M=input'whether or not use PID Controler with differential in advance:';for k=1:1:400timek=kts;%Linear modelyoutk=-den2y_1+num2u_5;kp=;kd=14;ki=;rink=signsin2pikts;rink=rink+sinpikts;errork=rink-youtk;ei=ei+errorkts;gama=;Td=kd/kp;c1=gamaTd/gamaTd+ts;c2=Td+ts/gamaTd+ts;c3=Td/gamaTd+ts;if M==1 %PID Control with differential in advanceudk=c1ud_1+c2youtk-c3y_1;uk=kperrork+udk+kiei;elseif M==2 %Simple PID Controluk=kperrork+kderrork-error_1/ts+kiei;endif uk>=110uk=110;endif uk<=-110uk=-110;end%Update parametersu_5=u_4;u_4=u_3;u_3=u_2;u_2=u_1;u_1=uk; y_3=y_2;y_2=y_1;y_1=youtk; error_2=error_1; error_1=errork; end figure1;plottime,rin,'r',time,yout,'b'; xlabel'times';ylabel'rin,yout'; figure2; plottime,u,'r'; xlabel'times';ylabel'u'; 其运行结果如表5所示;分析：微分先行PID 算法将微分算法放在前面,适用于给定量频繁升降的场合,可以避免给定值升降时引起的系统振荡,从而改善系统的动态特性;表5 例12仿真结果五、带死区的PID 控制算法某些系统为了避免控制作用过于频繁,消除由于频繁动作所引起的振荡,可采用带死区的PID 控制算法,控制算法为：⎩⎨⎧>≤=Bk e k e Bk e k e )()()(0)(,式中ek 为位置跟踪偏差,B 为可调的死区参数,具体可根据实际控制对象由试验确定;若B 太小,会使控制动作过于频繁,达不到稳定被控对象的目的；若B 太大,则系统将产生较大的滞后;Ex13 设被控对象为ss s s G 1047035.87523500)(23++=,采样周期为1ms,对象输出上有一个幅值为的正态分布的随机干扰信号;采用积分分离式PID 控制算法进行阶跃响应,取ξ=,死区参数B=,采用低通滤波器对对象输出信号进行滤波,滤波器为104.01)(+=s s Q ;仿真程序：;M=1时,采用一般积分分离式PID 控制算法,M=2时采用带死区的积分分离式PID控制算法;%PID Controler with dead zoneclear all;close all;ts=;sys=tf,1,,,0;dsys=c2dsys,ts,'z';num,den=tfdatadsys,'v'u_1=0;u_2=0;u_3=0;u_4=0;u_5=0;y_1=0;y_2=0;y_3=0;yy_1=0;error_1=0;error_2=0;ei=0;sys1=tf1,,1; %Low Freq Signal Filterdsys1=c2dsys1,ts,'tucsin';num1,den1=tfdatadsys1,'v';f_1=0;%M=1选择普通积分分离式PID,M=2选择带死区的积分分离式PID算法disp'M=1--Using common integration seperation PID Controler ,M=2-- Using integration seperation PID Controler with dead zone'M=input'whether or not use integration seperation PID Controler with dead zone:';for k=1:1:2000timek=kts;rink=1; %Step Signal%Linear modelyoutk=-den2y_1-den3y_2-den4y_3+num2u_1+...num3u_2+num4u_3Dk=rands1; %Disturbance signalyyoutk=youtk+Dk;%Low frequency filterfiltyk=-den12f_1+num11yyoutk+yy_1;errork=rink-filtyk;if abserrork<=ei=ei+errorkts;elseei=0;endkp=;ki=;kd=;uk=kperrork+kiei+kderrork-error_1/ts;if M==1uk=uk;elseif M==2 %Using Dead zoneif abserrork<=uk=0;endendif uk>=10uk=10;endif uk<=-10uk=-10;end%----------Return of PID parameters------------ rin_1=rink;u_3=u_2;u_2=u_1;u_1=uk;y_3=y_2;y_2=y_1;y_1=youtk;f_1=filtyk;yy_1=yyoutk;error_2=error_1;error_1=errork;endfigure1;subplot211;plottime,rin,'r',time,filty,'b';xlabel'times';ylabel'rin,yout';subplot212;plottime,u,'r';xlabel'times';ylabel'u';figure2;plottime,D,'r';xlabel'times';ylabel'Disturbance signal';其运行结果如表6所示;表6 例13仿真结果分析：在控制精度要求不太高,控制过程要求尽量平稳的场合,为了避免控制动作过于频繁,消除由此产生的震荡,可以认为设置一灵敏区;死区是一可调参数,参数太小,调节动作过于频繁,达不到稳定控制的目的,参数太大,又会产生很大的纯滞后;实验四纯滞后系统的控制算法一、纯滞后系统的Smith控制算法在工业过程控制中,许多被控对象具有纯滞后的性质;Smith提出了一中纯滞后补偿模型,其原理为,与PID控制器并接一个补偿环节,该补偿环节称为Smith预估器;被控对象传递函数为spesGτ-)(,模拟调节器的传递函数为)(sD,则系统的闭环传递函数为spspesGsDesGsDsττ--+=Φ)()(1)()()(;可见,闭环特征方程中出现了纯滞后环节,使系统稳定性降低,如果滞后时间τ足够大,系统将不稳定,这就是大延迟过程难于控制的本质;针对这一问题,Smith提出在调节器上反向并联一个预估模型,Smith预估控制系统等效图如图4-1所示;图4-1带Smith预估器的控制系统Ex14 设被控对象为160)(80+=-s e s G sp ,采用Smith 控制方法,在PI 控制中,取022.0,4==i p k k ,输入阶跃信号幅值为100;仿真程序：其运行结果如表1所示;分析：纯滞后性质常引起系统产生超调或者振荡,使系统的稳定性降低;s e τ-项在闭环控制回路之外,不影响系统的稳定性,仅将控制作用在时间坐标轴上推迟了一个时间τ,控制系统的过渡过程及其他性能指标都与对象特征为)(s G p 的完全相同,消除了纯滞后部分对控制系统的影响;表1 例14运行结果二、大林控制算法早在1968年,美国IBM 公司的大林就提出了一种不同常规PID 控制规律的新型算法, 即大林控制算法,该算法的最大特点是,将期望的闭环响应设计成一阶惯性加延迟,然后反过来得到能满足这种闭环响应的控制器;对于单回路控制系统,Dz 为数字控制器,Gz 为被控对象广义传递函数,则闭环系统传递函数为：)()(1)()()()()(z G z D z G z D z R z C z +==Φ则有)(1)()(1)()()(z z z G z E z U z D Φ-Φ==;如果能事先设定系统的闭环响应)(z Φ,则可得到控制器Dz;大林之处,通常的期望闭环响应是一阶惯性加纯延迟形式,其延迟时间等于对象的纯延迟时间τ：1)()()(+==Φ-s T e s R s C s sττ式中,τT 为闭环系统的时间常数,由此得到的控制率称为大林控制算法;Ex15 设被控对象为14.0)(76.0+=-s e s G sp ,采样时间为;期望的闭环响应设计为：115.0)()()(76.0+==Φ-s e s R s C s s仿真程序：;M=1时为采用大林控制算法,M=2时为采用普通PID算法;%Delay Control with Dalin Algorithmclear all;close all;ts=;%Plantsys1=tf1,,1,'inputdelay',;dsys1=c2dsys1,ts,'zoh';num1,den1=tfdatadsys1,'v';%Ideal closed loopsys2=tf1,,1,'inputdelay',;dsys2=c2dsys2,ts,'zoh';%Design Dalin controllerdsys=dsys2/dsys11-dsys2;num,den=tfdatadsys,'v';u_1=;u_2=;u_3=;u_4=;u_5=;y_1=;error_1=;error_2=;error_3=;ei=0;%选择是否使用大林算法,M=1使用,M=2不使用disp'M=1aUsing Dalin Method ,M=2-- Using common PID Controler' M=input'whether or not use Dalin Method:';for k=1:1:50timek=kts;rink=; %Tracing Step Signalyoutk=-den12y_1+num12u_2+num13u_3;errork=rink-youtk;if M==1 %Using Dalin Methoduk=num1errork+num2error_1+num3error_2+num4error_3-den3u_1-den4u_2-den5u_3-den6u_4-den7u_5/den2;elseif M==2 %Using PID Methodei=ei+errorkts;uk=errork+errork-error_1/ts+ei;end%----------Return of dalin parameters------------u_5=u_4;u_4=u_3;u_3=u_2;u_2=u_1;u_1=uk;y_1=youtk;error_3=error_2;error_2=error_1;error_1=errork;endplottime,rin,'b',time,yout,'r';xlabel'times';ylabel'rin,yout';其程序运行结果如表2所示;分析：大林控制算法,该算法的最大特点是,将期望的闭环响应设计成一阶惯性加延迟,然后反过来得到能满足这种闭环响应的控制器;大林算法很好的抑制了振铃现象;用大林算法设计的数字控制器系统输出超调量小甚至没有,并且允许有交长时间的调整时间;表2 大林算法与普通PID算法比较。

连续系统simulink状态空间建模分析方法程序设计连续系统simulink状态空间建模分析方法程序设计本课题基于对信号与信息处理课程中用matlab/simulink建模及应用分析滤波器问题的深入研究。

通过自身掌握的理论知识，主要以高阶连续系统（模拟滤波器）为例，并将其离散化，转化为离散系统，从而对离散系统处理。

用simulink状态空间函数模块建模，观察并分析波形。

其次，用matlab中的M文件编程，求解系统，绘制波形并进行频谱分析。

在本课题中，主要将连续系统转化为离散系统，再用计算机和matlab软件进行研究，用simulink对高阶离散系统建模，并设置模块参数，自定义函数为正弦波的叠加，传输函数的相关参数后运行并进行频谱分析，使信号的性态都能得到处理和研究。

通过编程，求解高阶离散系统的零输入响应，零状态响应和完全响应，求解实际生活中的各种问题，改变参数并对信号进行适当的频谱分析。

目录引言 (2)1 状态空间分析方法的概述 (2)2 快速创建LTI连续系统状态空间模型的方法 (3)2.1 创建LTI连续系统传递函数的方法 (3)2.2 构造描述LTI连续系统的状态空间模型矩阵 (4)3 利用simulink状态空间建模求解LTI系统数值解的思路 (5)3.1 MATLAB编程设计并描述低通数字滤波器 (5)3.2 创建系统的simulink状态空间模型 (7)3.3 模块内部参数设置及数据存储 (8)4 利用simulink状态空间建模求解LTI系统的优缺点 (9)5 连续系统simulink状态空间建模分析方法程序设计的思路 (10)5.1 调用模型文件及编程求解系统响应 (10)5.2 分析系统的频谱与相位 (11)6 状态空间分析方法的应用实例 (13)6.1 分析求解低阶电路系统 (13)6.2 设计分析滤波器系统...................................... 14 7 结束语 ..................................................... 14 引言随着科学技术的发展，系统的组成也日益复杂。

M 文件和Simulink 求解连续微分系统实例分析Matlab 自带的一个S 函数源代码： D:\MATLAB7\toolbox\simulink\blocks\sfuntmpl例1. 常微分方程（Lorenze 混沌系统）：112322331223,,,x b x x x x a x a x x x x c x x=-+⎧⎪=-+⎨⎪=-+-⎩ （1） 其中10,28,8/3a r b ===。

（1） m 文件实现：文件名为exam1.mfunction exam1x0=[0;0;1e-3];[t,x]=ode45(@lorenzfun,[0,100],x0); figure(1) plot(t,x) figure(2)plot3(x(:,1),x(:,2),x(:,3))%----------------------------------function dx=lorenzfun(t,x)a=10;c=28;b=8/3; dx=zeros(3,1);dx(1)=-b*x(1)+x(2)*x(3); dx(2)=-a*x(2)+10*x(3);dx(3)=-x(1)*x(2)+c*x(2)-x(3);(2) (I)Simulink模块实现：（见lorenzblok）其中三个积分模块的初始值设置与exam1相同，仿真时长为100s。

精度设置：Simulation--Configuration Parameters—Relative tolerance, 1e-3改为1e-5（试试不作此修改的结果比较）。

运行后双击示波器scope后可看到：或在matlab命令窗口输入画图命令：>> plot(tout,yout) >>plot3(yout(:,2),yout(:,3),yout(:,1))(II)Simulink 向量模块实现：(见lorenzevector)画图语句：>>plot(t,x) >>plot3(x(:,1),x(:,2),x(:,3))(3)Simulink 中S 函数的实现：（见lorenzsfun 和lorenzsystem.m ）例2. 常时滞微分方程:'11'212'32(1),(1)(0.2),(),y y t y y t y t y y t =-=-+-= （2）（1） m 文件需调用dde23来求解：（见exam2.m ）function exam2sol = dde23('exam1f',[1, 0.2],ones(3,1),[0, 5]);plot(sol.x,sol.y); title('Example 2') xlabel('time t'); ylabel('y(t)');%----------------------------------------------- function v = exam1f(t,y,Z) ylag1 = Z(:,1); ylag2 = Z(:,2); v = zeros(3,1);v(1) = ylag1(1);v(2) = ylag1(1) + ylag2(2); v(3) = y(2);Example 2time ty (t )(2)Simulink 中S 函数来实现：（见exam2sfun 和exam2mfun.m ） 注：用Simulink 中S 函数求解时滞微分方程的核心思想在于：将时滞变量作为S 函数的外部输入。