1、2018年9月26日星期三,1,Matlab 2008b下载地址：ed2k:/|file|%E7%9F%A9%E9%98%B5%E5%AE%9E%E9%AA%8C%E5%AE%A4.Mathworks.Matlab.R2008b.DVD.ISO-TBE.iso|4306182144|A2EDE3392561D89CE8EA3C37CFE2429F|h=HREWLFRSPC6YM65XD52FMVORHGGEZHVL|/网页：http:/ 考虑带有一阶延迟的多变量传递函数矩阵,用MATLAB命令进行阶跃响应,c5eg3a,step响应是在两路输入单独作用下分别得到的，表明Gij(s) (ij)对
2、输入有很强的响应，说明系统有较强的耦合。为了对两路信号分别设计控制器，希望减小这种耦合。为此考虑对系统进行补偿。,2018年9月26日星期三,3, c5eg3b,例3 考虑带有一阶延迟的多变量传递函数矩阵,考虑补偿矩阵Kp,并对G(s)*Kp进行单位阶跃响应。由于对象中含有时滞，不能利用矩阵乘法，对时滞要用pade逼近，或者在simulink下仿真。,2018年9月26日星期三,4,考虑用simulink仿真，并与pade逼近的结果进行比较,1.建立模型,2.修改模型参数,2.1两个阶跃信号起始时间改为0，终值时间分别为u1,u2,2.2 修改矩阵增益为 0.1134, 0.924; 0.33
3、78, -0.318,2.3 修改传递函数,2018年9月26日星期三,5,考虑用simulink仿真，并与pade逼近的结果进行比较,1.建立模型,2.修改模型参数,2.4 将时滞分别改为 0.72, 0.30, 1.29.,3.用simulink模型仿真,2018年9月26日星期三,6,考虑用simulink仿真，并与pade逼近的结果进行比较,1.建立模型,2.修改模型参数,3.用simulink模型仿真, c5eg3c,带有Pade逼近的第一输入通道输出,带有Pade逼近的第二输入通道输出,Simulink仿真的第一输入通道输出,Simulink仿真的第二输入通道输出,2018年9月2
4、6日星期三,7,例4 计算机控制系统仿真,考虑如图所示的计算机控制系统，其中控制器模型D(z)是离散模型，采样周期为Ts，ZOH为零阶保持器，而受控对象模型G(s)为连续模型:,对这样的系统直接写成微分方程的形式再进行仿真的方法是行不通的，因为其中既有连续环节，又有离散环节。所以用simulink仿真。,模型中参数a, T 需要用户确定，,2018年9月26日星期三,8,Simulink仿真模型如下：,模型中参数 a, T 需要用户确定，而 z1, p1, K 需要由控制器模型计算。,stairs(x,y)画出由向量 x, y 确定的阶梯图形,2018年9月26日星期三,9,仿真结果如下,T
5、= 0.2s 时,T = 1s 时,2018年9月26日星期三,10,对于该例，也可采用离散与连续传递函数转换的方法，在给定采样周期下获得离散传递函数。, c5eg4b,Zero/pole/gain: 0.018187 (z+0.9934) (z-0.9802)-(z-0.9802) (z2 - 1.801z + 0.8368) Sampling time: 0.2,2018年9月26日星期三,11,例5 时变系统的仿真,考虑如图所示控制系统模型，其中控制器参数Kp=200, Ki=10, 饱和非线性宽度为 2，受控对象为时变模型，由以下微分方程给出,试分析该系统的阶跃响应曲线。,设,，则可将
6、系统写成,2018年9月26日星期三,12,系统的 simulink 模型如下,注意：1. 函数模块中自变量是用输入信号u表示的，u是由时钟模块产生； 2. 生成时变部分模型与状态变量用乘法器相成即可。,2018年9月26日星期三,13,c5eg5a,2018年9月26日星期三,14,例6 多采样速率系统的仿真,如图所示双环电机控制系统，内环为电流环，采样周期为T10.001s，控制器模型为D1(z)=(0.0967z0.0965)/(z1), 控制器外环采样周期为T2 = 0.01s，控制器模型为D2(z)=(5.2812z5.2725)/(z1)。,搭建simulink仿真模型，并分析阶跃
7、响应。,2018年9月26日星期三,15,Simulingk模型为, t,x,y = sim(c5eg6, 2); plot(t,y),2018年9月26日星期三,16,例7. 系统的脉冲响应分析,在simulink里并没有提供单位脉冲信号模块，所以可以用阶跃模块来近似，如令阶跃初始值为1/a，阶跃时间(即阶跃宽度)为a，阶跃终止时为0, 如图。,阶跃初值,阶跃宽度,2018年9月26日星期三,17,例如将如图所示系统用近似单位脉冲作为输入，观察动态响应,2018年9月26日星期三,18,指定输入向量: -3,-2,-1, 2,3,4 ,分段线性的静态非线性环节,指定输出向量: -1,-1,0
8、,0,1,1,-3 -2 -1 0 1 2 3 4,1 01,2018年9月26日星期三,19,斜坡输入响应,分段线性函数响应,2018年9月26日星期三,20,分解,例8. 用simulink搭建如图所示多值非线性环节,当u增加时，取y1(u)；当u减小时，取y2(u)。Simulink模型如图。,修改模块参数,2018年9月26日星期三,21,例8. 用simulink搭建如图所示多值非线性环节,分解,记忆模块，保存上一步的值uk1 .初值为0,将第1端口输入值 ( 当前值uk )与第2端口输入值 ( 上一步的值 uk1) 比较, uk uk1时输出1，否则为0。,当第2输入为1时(大于设
9、定的阈值0.5)，模块输出取第1输入，当第2输入为0时(小于设定的阈值0.5) ，模块输出去第3输入。,2018年9月26日星期三,22,例8. 用simulink搭建如图所示多值非线性环节,对正弦输入u(t)= asin t的响应,a=2; t, x1, y1=sim(c5eg8, 10);a=4; t, x2, y2=sim(c5eg8, 10);a=8; t, x3, y3=sim(c5eg8, 10);plot(t, y1, y2, y3),2018年9月26日星期三,23,例9. 类似可以建立如下多值非线性,x1=-3,-2,-1,2,3,4, y1=-1,-1,0,0,1,1x2=
10、-3,-2,-1,1,2,3, y2=-1,-1,0,0,1,1,对正弦输入u(t)= asin t的响应,a=2; t, x1, y1=sim(c5eg8, 10);a=4; t, x2, y2=sim(c5eg8, 10);a=8; t, x3, y3=sim(c5eg8, 10);plot(t, y1, y2, y3),2018年9月26日星期三,24,非线性系统极限环的仿真研究,例10. 考查如图所示的典型非线性模型，其中非线性环节为例8中的回环非线性。,用simulink搭建仿真模型如图。积分器模块的初值设置为1，就可获得系统的自激振荡。,2018年9月26日星期三,25, opt=
11、simset(RelTol,1e-8); t,x,y=sim(c5eg10,40,opt); subplot(121),plot(t,y) subplot(122),plot(y(:,1),y(:,2),2018年9月26日星期三,26,十二. 编写M语言S-函数,S-函数 S-函数是Simulink模块的计算机语言描述，它可以用MATLAB、C、C+、Ada或Fortran语言编写。MATLAB 语言的 S-函数可编译为M-文件。,S-函数使用特定的调用语法与Simulink中的方程求解器相互作用，接收求解器中的信息，并对求解器的指令作出响应。,S-函数的格式非常通用，可用在连续、离散、以及
12、混杂系统中。 用户可以用S-函数实现自己的算法，并可形成Simulink模块。,在模型中利用S-函数 1).首先将 S-函数模块从Simulink /User-Defined Functions模型库中拖拽到模型窗口;,2).打开模型对话框, 在“S-Function name:” 栏内填入S-函数名称，在“S-Function parameters:”栏内填入S-函数的参数。,2018年9月26日星期三,27,function sys,x0,str,ts=mysfun(t,x,flag, p1,p1)% mysfun M-file% switch (flag).,S-函数名称，填入对话框S-
13、function name栏,p1, p2: S-函数参数，在S-function parameters栏中赋值,输入S-函数名,允许用户向S-函数传递参数, 用户必须知道参数的顺序、类型, 它们可以是常数或变量, 不同参数值之间用逗号隔开,2018年9月26日星期三,28,S-函数的工作方式Simulnk模块的数学含义,其中，x=xc+ xd , Ts为离散状态的采样时间。,仿真过程初始化阶段： Simulink将库模块结合到模型中, 传递信号宽度, 确定数据类型, 采样时间, 计算模块参数, 确定模块执行顺序并分配内存。,仿真循环阶段： 每通过一个仿真循环称为一个仿真步, 在每一仿真步内,
14、 Simulink 按照初始化确定的顺序执行模型的每一模块, 对于每一模块, Simulink 调用相应函数计算状态, 导数, 并计算当前采样时间的输出直到仿真结束。,2018年9月26日星期三,29,S-函数回调方法 S-函数由一组S-函数回调方法组成, 这些方法在每个仿真阶段执行不同任务, 这些任务包括: 初始化 -初始化SimStruct, 这是一个包含S- 函数信息的仿真结构; -设置输入和输出端口的个数和维 数; -设置模块采样时间; -分配空间和size数组。 计算下一个采样时间点(适于变采 样时间模块)。 以最大时间步骤计算输出。 以最大时间步更新离散状态。 积分(适于连续状态)
15、。,仿真步骤示意图,2018年9月26日星期三,30,M文件S-函数的实现 M文件S-函数包含下列函数:,sys, x0,str,ts=f(x,t,u,flag,p1,p2,),f: S-函数的名称; x: S-函数对应模块的状态向量; t: 当前时间; u: 模块输入; flag: 标识要执行的任务(子函数);p1,p2: 模块参数;,sys: 其含义根据标志flag 取值不 同返回不同子函数的返回值; x0: 初始化状态向量; str: 空矩阵(留做将来使用); ts: 采样时间.,S-函数的内容,2018年9月26日星期三,31,S-函数的一些概念,直接馈通:模块的输出(或变采样时间的采
16、样时间)直接接收输入端口值的控制。如 输出(flag=3)是输入 u 的函数; 变采样时间S-函数中的下一次采样 时间点(flag=4)的确定要访问输入u.,动态设置数组: M文件 S-函数只有一个输入端口, 可以接收一个向量信号, 但信号宽度是可以不同的, 是可以动态指定的 (可以在sizes结构中指定相应的属性值为1 ), 也可以由输入的维数确定输入端口信号宽度, 连续和/或离散状态的数目. 如果使用length(u)函,数调用S-函数, 则用户可以确定模块实际输入的宽度, 如指定宽度为零, 则输入端口被从S-函数模块中删除. 例如:图中有两个相同的S-函数模块, 上面的模块由一个带有三个
17、输出分量的Mux向量模块驱动, 下面的S-函数模块由一个带有标量输出的模块驱动.,2018年9月26日星期三,32,采样时间和偏移量的设置 在M文件S-函数中, 指定采样时间给出了如下选项., 离散采样时间适于用户的S-函数模块的行为带有离散时间间隔, 用户可以定义采样时间以便控制simulink何时调用S-函数模块, 也可以定义延迟每个采样点的偏移量, 偏移量的数值不超过对应的采样时间. 采样点时刻由下式确定, 连续采样时间适于连续状态和/或非采样过零的S-函数, 输出以最小时间步改变., 连续但在最小时间步固定的采样时间适于在每一最大仿真步上执行, 但在最小时间步上不改变数值的S-函数.,
18、下一采样时刻=(n*采样间隔)+偏移量n 从0开始, 表示当前仿真步数. 在每个采样时刻Simulink调用S-函数的mdlOutput函数和mdlUpdate函数来计算输出和更新状态., 变采样时间两次采样之间的时间间隔是可变的, 在每次仿真步开始时需计算下一仿真时刻., 继承采样时间用户可以指定模块的采样时间为inherited(继承性). 对于一个模块, 它可以继承下列模块的采样时间: -驱动模块; -目标模块; -系统中最快的采样时间.相应地, 在M文件S-函数中设置采样时间为1.,2018年9月26日星期三,33,采样时间和偏移量的指定形式 ts = 采样时间1, 偏移量1; 采样时
19、间2, 偏移量2; 采样时间k, 偏移量k S-函数可以使用单速率系统和多速率系统, 多速率系统有多个采样时间.,ts=-1.0, Offset, 从驱动模块继承采 样时间, 但在最小时间步内不随 输入改变.,有效的采样时间对为:ts=0.0, 0.0, 连续采样时间;ts=0.0, 1.0, 连续但在最小时间步 固定的采样时间;ts=-2.0, 0.0, 变采样时间;ts=Period, Offset, 离散采样时间;ts=-1.0, 0.0, 从驱动模块继承采样 时间;,对于多速率S-函数, 每个任务以不同的速率执行, 这时, ts应该以采样时间上升的顺序指定用户 S-函数中使用的采样速率
20、. 如任务i 的采样时间为: i , i= 1, 2, 且 1 2 2 1ts = 1, ofset1; 2, ofset2则Simulink在下列时刻执行S-函数,0, ofset1, ofset2, ofset1 +1, ofset2 +2 , ofset1 +21, ofset2 +22, 且在ofset1 +k 1时执行任务1, ofset2 +k 2时执行任务2.,2018年9月26日星期三,34,M文件S-函数模板 为使Simulink识别M文件S-函数, 用户必须提供S-函数的某些特定信息, 在初始化阶段(flag=0), S-函数调用子函数mdlInitializeSizes，
21、对一个sizes结构变量初始化, sizes 结构包含了S-函数的特定信息.,Sizes的结构属性,2018年9月26日星期三,35,S-函数模板,M-文件 S函数模板sfuntmpl.m 位于 matlabroot/toolbox/simulink/blocks,我们可以通过修改该模板来建立我们自己的S-函数,双击可以打开模板,2018年9月26日星期三,36,function sys,x0,str,ts = sfuntmpl(t,x,u,flag)% The following outlines the general structure of an S-function.switch f
22、lag,case 0, sys,x0,str,ts=mdlInitializeSizes; % Initialization case 1, sys=mdlDerivatives(t,x,u); % Derivativescase 2, sys=mdlUpdate(t,x,u); % Update case 3, sys=mdlOutputs(t,x,u); % Outputscase 4, sys=mdlGetTimeOfNextVarHit(t,x,u); % Get Time Of Next VarHitcase 9, sys=mdlTerminate(t,x,u); % Termina
23、teotherwise error(Unhandled flag = ,num2str(flag); % Unexpected flagsEnd%=,2018年9月26日星期三,37,%=% mdlInitializeSizes% Return the sizes, initial conditions, and sample times for the S-function.%=function sys,x0,str,ts=mdlInitializeSizessizes = simsizes;sizes.NumContStates = 0;sizes.NumDiscStates = 0;si
24、zes.NumOutputs = 0;sizes.NumInputs = 0;sizes.DirFeedthrough = 1;sizes.NumSampleTimes = 1; % at least one sample time is neededsys = simsizes(sizes);x0 = ; % initialize the initial conditionsstr = ; % str is always an empty matrixts = 0 0; % initialize the array of sample times%=,2018年9月26日星期三,38,%=%
25、 mdlDerivatives% Return the derivatives for the continuous states.%=%function sys=mdlDerivatives(t,x,u)sys = ; %=% mdlUpdate% Handle discrete state updates, sample time hits, and major time step% requirements.%=%function sys=mdlUpdate(t,x,u)sys = ; %=,2018年9月26日星期三,39,%=% mdlOutputs% Return the bloc
26、k outputs.%=%function sys=mdlOutputs(t,x,u)sys = ; %=% mdlGetTimeOfNextVarHit% Return the time of the next hit for this block. Note that the result is% absolute time. Note that this function is only used when you specify a% variable discrete-time sample time -2 0 in the sample time array in% mdlInit
27、ializeSizes.%=%function sys=mdlGetTimeOfNextVarHit(t,x,u)sampleTime = 1; % Example, set the next hit to be one second later.sys = t + sampleTime;%=,2018年9月26日星期三,40,%=% mdlTerminate% Perform any end of simulation tasks.%=%function sys=mdlTerminate(t,x,u)sys = ;% end mdlTerminate,2018年9月26日星期三,41,静态系
28、统S-函数,例1 y = 2u,function sys,x0,str,ts = timestwo(t,x,u,flag)%TIMESTWO S-function whose output is % two times its input.switch flag, case 0 sys,x0,str,ts=mdlInitializeSizes; case 3 sys=mdlOutputs(t,x,u); case 1, 2, 4, 9 sys=; otherwise error(Unhandled flag = ,num2str(flag);end%=,function sys,x0,str,
29、ts = mdlInitializeSizes()sizes = simsizes;sizes.NumContStates = 0;sizes.NumDiscStates = 0;sizes.NumOutputs = -1; sizes.NumInputs = -1; sizes.DirFeedthrough = 1; sizes.NumSampleTimes = 1;sys = simsizes(sizes);str = ;x0 = ;ts = -1 0; % inherited sample time%=function sys = mdlOutputs(t,x,u)sys = u * 2
30、;,timetwo.m S-函数代码如下,2018年9月26日星期三,42,进行了封装,输出幅值为1:2的正弦波,输出为幅值为1:3, 周期为1, 宽度为0.5s的矩形波,2018年9月26日星期三,43,例2 y = ku , k为用户指定的静态增益,sfun_vargain.m S-函数代码如下,function sys,x0,str,ts = sfun_vargain(t,x,u,flag,gain)switch flag, case 0 sys,x0,str,ts=mdlInitializeSizes; case 3 sys=mdlOutputs(t,x,u,gain); case 1
31、, 2, 4, 9 sys=; otherwise error(Unhandled flag = ,num2str(flag);end%=function sys,x0,str,ts = mdlInitializeSizes()sizes = simsizes;,sizes.NumContStates = 0;sizes.NumDiscStates = 0;sizes.NumOutputs= -1; % 动态指定sizes.NumInputs= -1; % 动态指定sizes.DirFeedthrough = 1; % 有直接馈通sizes.NumSampleTimes = 1;sys = s
32、imsizes(sizes);str = ;x0 = ;ts = -1 0; % 继承采样时间%=function sys = mdlOutputs(t,x,u,gain)sys = u * gain;,2018年9月26日星期三,44,打开对话框, 选择增益,填入所选增益, 如4,仿真,2018年9月26日星期三,45,连续系统S-函数,例3. 建立一个积分器, 它的初值作为用户输入. 状态方程为,sfun_int.m S-函数代码如下,function sys,x0,str,ts = sfun_int(t,x,u,flag,initial_state)switch flag, case 0
33、, sys,x0,str,ts= mdlInitializeSizes(initial_state); case 1, sys=mdlDerivatives(t,x,u); case 3, sys=mdlOutputs(t,x,u); case 2,4,9 sys=;otherwise error(Unhandled flag = ,num2str(flag);end,%=function sys,x0,str,ts= mdlInitializeSizes(initial_state)sizes = simsizes;sizes.NumContStates = 1;sizes.NumDiscS
34、tates = 0;sizes.NumOutputs = 1;sizes.NumInputs = 1;sizes.DirFeedthrough = 0;sizes.NumSampleTimes = 1; sys = simsizes(sizes);x0 = initial_state;str = ;ts = 0 0;%=function sys=mdlDerivatives(t,x,u)sys = u;%=function sys=mdlOutputs(t,x,u)sys = x;,2018年9月26日星期三,46,打开对话框, 选择初始状态,填入所选初始状态, 如5,仿真,2018年9月26
35、日星期三,47,例4. 蹦极跳系统的微分方程可写为:,k: 弹性系数 m: 物体质量 g: 重力加速度a1, a2: 空气阻尼系数,令 , h为距离地面的距离, 则系统方程可写为:,构造S-函数, 取k =5; a1=1; a2=1; g=10. 输入长度(x1的初始值), 质量m和距离地面的距离h.,2018年9月26日星期三,48,sfun_bungee.m S-函数代码如下,function sys,x0,str,ts = sfun_bungee(t,x,u,flag,len,weight,dist_ground)% bungee jumping sfunction:k=5;a1=1;a
36、2=1;g=10;% length,weight,distance from ground are input parametersk=5;a1=1;a2=1;g=10;switch flag, case 0, sys,x0,str,ts=mdlInitializeSizes(len); case 1, sys=mdlDerivatives(t,x,u,weight,k,a1,a2,g); case 2, sys=mdlUpdate(t,x,u); case 3, sys=mdlOutputs(t,x,u,dist_ground); case 4, sys=mdlGetTimeOfNextVa
37、rHit(t,x,u); case 9, sys=mdlTerminate(t,x,u); otherwise error(Unhandled flag = ,num2str(flag);end,输入参数,输入参数,输入参数,设定参数,设定参数,输入参数,2018年9月26日星期三,49,%=function sys,x0,str,ts=mdlInitializeSizes(len)sizes = simsizes;sizes.NumContStates = 2;sizes.NumDiscStates = 0;sizes.NumOutputs = 1;sizes.NumInputs = 0;s
38、izes.DirFeedthrough = 0;sizes.NumSampleTimes = 1; % at least one sample time is neededsys = simsizes(sizes);x0 = -len;0;str = ;ts = 0 0;%=,2018年9月26日星期三,50,%=function sys=mdlDerivatives(t,x,u,weight,k,a1,a2,g)if x(1)0 b=0;else b=-k*x(1);endx1dot=x(2);x2dot=1/weight*(weight*g+b-a1*x(2)-a2*abs(x(2)*x(
39、2);sys = x1dot;x2dot;% end mdlDerivatives%=function sys=mdlUpdate(t,x,u)sys = ;% end mdlUpdate%=function sys=mdlOutputs(t,x,u,dist_ground)sys = dist_ground-x(1);% end mdlOutputs%=,2018年9月26日星期三,51,%=function sys=mdlGetTimeOfNextVarHit(t,x,u)sampleTime = 1; % Example, set the next hit to be one second later.sys = t + sampleTime;% end mdlGetTimeOfNextVarHit%=function sys=mdlTerminate(t,x,u)sys = ;% end mdlTerminate,
