非连续系统的Simulink仿真方法研究

1、第18卷第7期 系统 仿 真 学 报© V ol. 18 No. 72006年7月 Journal of System Simulation July, 2006非连续系统的Simulink 仿真方法研究向 博1,高丙团1,张晓华1,胡广大2(1.哈尔滨工业大学电气工程学院,哈尔滨 150001;2.哈尔滨工业大学航天学院,哈尔滨 150001 摘 要:针对滑模变结构控制、最优控制等控制系统中出现的由于不连续性问题而直接影响Simulink 仿真速度和仿真效果的情况,分析了该现象产生的原因,并提出了相应的解决策略。结合吊车防摆变结构控制系统,采用提出的解决方案进行了Simulink

2、仿真实验研究,仿真实验的结果证明所提出的策略对于解决这类不连续性仿真问题具有较好的效果。最后,给出了解决和避免这种问题的几点建议。 关键词:非连续系统;Simulink 仿真;滑模变结构控制;吊车中图分类号:TP391.9 文献标识码:A 文章编号:1004-731X (2006 07-1750-05Simulink Simulation Research of Uncontinuous SystemXIANG Bo, GAO Bing-tuan, ZHANG Xiao-hua, HU Guang-da(1.School of Electrical engineering, Harbin In

3、stitute of Technology, Harbin 150001, China;2.School of Astronautics, Harbin Institute of Technology, Harbin 150001, ChinaAbstract: The Simulink simulation speed and effect are decreased because of the uncontinuous functions in the uncontinuous systems, such as sliding mode variable structure contro

4、l systems, optimal control systems and so on. The reasons caused the problem were analyzed in detail and corresponding resolving methods were also developed. Based on the presented methods, the Simulink simulation on variable structure control of an overhead crane system was studied, and the simulat

5、ion results demonstrate the effectiveness of the prompted Simulink simulation schemes. Finally, some suggestions on the problem were summed up.Key words: uncontinuous system; Simulink simulation; sliding mode variable structure control; overhead crane引 言 Simulink 是MathWorks 公司于1990年推出的产品,是用于MATLAB 下

6、建立系统框图和仿真的环境。该环境刚推出时的名字叫Simulab ,由于其名字很类似于当时的一个很著名的语言Simula 语言,所以次年更名为Simulink 。从名字上看,立即就能看出该程序有两层含义1,首先,“Simu ”一词表示它可以用于计算机仿真,而“Link ”一词表明它能进行系统连接,即把一系列模块连接起来,构成复杂的系统模型。正是由于它的这两大功能和特色,使得它成为仿真领域首选的计算机环境。Simulink 是一个用来对动态系统进行建模、仿真分析的软件包,它支持连续、离散及两者混合的线性和非线性系统,也支持具有多种采样速率的多速率系统。其利用数学模型来获取系统的一些重要特性参数,这

7、些数学模型通常是以时间为变量的常微分方程来描述,并用数值方法进行计算机仿真求解。所以,Simulink 软件包常被用于运动学仿真来反复求解 机构运动约束方程2,通过积分获得最终的速度(或加速度,收稿日期:2005-05-25 修回日期:2006-04-03作者简介:向博(1980-,男,四川广元人,硕士生,研究方向为基于滑模变结构控制的吊车防摆系统;张晓华(1961-,男,黑龙江哈尔滨人,教授,博导,哈尔滨工业大学电气工程学院常务副院长,研究方向为管道机器人,MATLAB 仿真技术。从而确定机构运动的位置(和速度。图1为运动学仿真的流程图,它描述了诸如MATLAB 仿真工具箱Simulink

8、仿真软件包的仿真步骤和流程2。图1 运动学仿真流程框图利用Simulink 仿真软件包来求解机构运动,比通常求解机构运动的方法具有许多优点。最显著的优点是位置求解问题可以用软件包隐式求解,使用者只需为仿真提供适当的初始条件即可以求解机构在任意时刻的位置问题。但是,如果在仿真模型中一旦出现了非连续模块,则可能出现系统仿真停滞的情形。本文针对这一现象,以“基于滑模变结构控制的桥式吊车防摆系统”为例,详细分析了该现象产生的原因,并提出了相应的解决策略。1 问题的提出本文以“基于滑模变结构控制的吊车防摆系统”为例,来研究非连续系统的Simulink 仿真方法。2006年7月 向 博,等:非连续系统的S

9、imulink 仿真方法研究 July, 2006根据文献3,吊车的物理模型如图2所示:x 图2 吊车的物理模型重物通过绳索与小车相连,小车在行走电机的水平驱动力f 的作用下在水平轨道上运动。小车的质量为M ,重物的质量为m ,绳索的长度为l ,小车与水平轨道的摩擦阻尼系数为D ,g 为重力加速度时间常数。利用拉格朗日方程,我们可以得到吊车的数学模型为: f Mx Dx mg x l g =+=+(1 1.1 滑模变结构控制器的设计设1x =,2x = ,1l m =,29.8/g m s =,并取控制:a x =则根据式(1,可写出如下二阶系统:12219.8xx x x a =+ (2 取

10、滑模面函数为:112s c x x =+ 10c > (3采用等速趋近率控制策略:sgn(dss dt= 0> (4 设李亚甫诺夫函数:212V s =(5 则可得:sgn(0Vs s s s s =< (6 由于在切换面邻域内函数式(5是正定义的,而式(6中2s 的导数是负半定义的。所以,在s =0附近V 是一个非增函数,因此,如果满足条件式(6,则定义函数式(5是系统的一个条件李亚普诺夫函数。系统本身也就稳定于条件0s =。系统将沿滑模面在t 趋于无穷大时归于平衡点。由(2、(3、(4式变换可得:1121129.8sgn(a x c x c x x =+(7 式(2、(7

11、即构成了具有自动防摆功能的二维滑模变结构控制器。1.2 系统仿真针对前面所设计的滑模变结构控制系统进行仿真,搭建的仿真模型如图3所示。在此模型中,我们全部采用Simulink 工具箱中的模块搭建模型。图3 滑模变结构控制器仿真模型取11c =,1=,1(01x =,2(00x =,仿真时间T=10s , 采用Variable-step 的ODE45算法。仿真结果如图4所示:图4 系统时域仿真结果由上图可以看出,当系统仿到1s 左右时,速度非常慢,停滞不前,出现了图形未完全仿真完成的情形,这显然是我们不希望看到的。2 问题的分析为了更好的分析该问题产生的原因以及提出相应的解决策略,首先,我们认为

12、有必要从Simulink 仿真运行原理开始进行分析。2.1 Simulink 仿真运行原理对Simulink 系统模型的仿真主要包含两个阶段4:初始化阶段和运行阶段1.初始化阶段在初始化阶段主要要完成以下工作:(1 每个模块的所有参数都传递给MATLAB 进行求值,得到的数值作为实际的参数使用;(2 展开模型的层次结构,每个子系统被他们所包含的模块替代,带有触发和使能模块的子系统被视为原子单元进行处理;2006年7月系统仿真学报 July, 2006(3 检查信号的宽度和模块的连接情况,提取状态和输入/输出依赖关系方面的信息,确定模块的更新顺序;(4 确定状态的初值和采样时间;完成这些工作后就

13、可以进行仿真了。2.运行阶段初始化之后,仿真进入运行阶段。仿真是由求解器控制的,它计算模块的输出,更新离散状态,计算连续状态。在采用变步长求解器时,求解器还会确定时间步长。计算连续状态包含下面两个步骤:(1 首先,求解器为待更新的系统提供当前状态、时间和输入值。反过来,求解器需要状态导数的值;(2 然后,求解器对状态的导数进行积分,计算新的状态的值。状态计算完成后,再进行一次模块的输出更新,这时,一些模块可能会发出过零的警告,促使求解器探测出发生过零的准确时间。Simulink的仿真过程是在Simulink求解器和系统相互作用之下完成的,系统和求解器在仿真过程中的对话作用如图5所示。 图5 系

14、统与求解器之间的对话在图5中,求解器的作用是传递模块的输出,对状态导数进行积分,并确定采样时间。系统的作用是计算模块的输出,状态进行更新,计算状态的导数,生成过零事件、从求解器传递给系统的信息包括时间、输入和当前状态;反过来,系统为求解器提供模块的输出、状态的更新和状态的导数。2.2 过零检测过零检测5是通过在系统和求解器之间建立对话的方式工作,对话包含的一个内容是事件通知,即系统告知求解器在前一时间步发生了个事件。事件由过零表示,过零在下列两个条件下产生:(1 信号在上一个时间步改变了符号(含变为O和离开O;(2 模块在上一个时间步改变了模式(例如积分器进入了饱和区段。过零是一个重要的事件,

15、表征系统中的不连续性。例如,响应中的跳变。如果仿真中不对过零进行检测,可能会导致不准确的仿真结果。当采用变步长求解器时,Simulink 能够检测到过零(使用固定步长的求解器,Simulink不检测过零。当一个模块通知系统前一时间步发生了过零,变步长求解器就会缩小步长,即便绝对误差和相对误差是可接受的。缩小步长的目的是判定事件发生的准确时间。当然,这样会降低仿真的速度,但这样做对有些模块来讲是至关重要和必要的,因为这些模块的输出可能表示了一个物理值,它的零值有着重要的意义。事实上,只有少量的模块能够发出过零事件通知。每个模块发出专属于自己的事件通知,而且可能与不止一个类型的事件发生关联。能够产

16、生过零通知的Simulink模块如图6所示:图6 能够产生过零通知的SIMULINK模块各个模块的过零的类型是有差异的。例如,abs模块在输入改变符号时产生个事件,而saturation模块能够生成两个不同的过零,一个用于下饱和,一个用于上饱和。signals&system库中的HitCrossing模块输入穿过零点时产生一个过零,可以用来为不带过零能力的模块提供过零检测的能力。一些过零只是用来通知求解器模式已经发生了改变,另外一些则与信号相关,用于触发其他模块。触发包含三种类型:(1 上升沿信号上升到或穿过零或者信号离开零变正;(2 下降沿信号下降到或穿过零或者信号离开零变负;(3

17、双边沿上升或下降两者之一发生。注意,Fcn模块不支持过零,结果一些拐角点被漏掉了。2.3 算法的选择Simulink提供了许多微分方程的解法6,如Runge-Kutta 法、Adams法、Gear法、Euler法、Linsim法等等。需要说明的是,对于不同的仿真问题,需要按照系统的特性,选择不同的仿真方法。一般来说:Runge-Kutta方法(包含rk45、rk23变步长方法适合于高度非线性或不连续的系统,不适合于刚性系统(即有快变特性又有慢变特性的系统;Adams方法适合于非线性小、时间常数变化小的系统;Gear方法是专门用于刚性系统仿真的方法,但对非刚性系统较差;Euler方法比较差,尽量

18、避免使用;Linsim方法适合于近似线性的系统,对线性刚性系统有很大的优越性。一般来说,使用变步长的自适应算法是较好的选择。这类算法会依照给定的精确度在各积分区内自适应地寻找各自的最大步长进行积分,从而使得速率最高。Simulink的变步2006年7月向博,等:非连续系统的Simulink仿真方法研究July, 2006长解法能够把积分段分到足够细以至到满足精度要求的解。2.4 系统Simulink仿真停滞原因分析综上所述,我们可以发现:由于系统存在不连续模块signum,所以,当系统于1s左右到达滑模面(s=0时,signum 模块向系统发出过零通知。而当采用变步长求解器时, Simulin

19、k能够检测到过零。所以,当Simulink检测到过零时,便自动缩小步长,可是,在下一仿真步里,系统继续过零。如图7和如图8所示,因为滑模面在1s处不能正常归零,所以signum模块就反复过零,同时一直向求解器发出过零通知。求解器便相应的一直不停的缩小步长。 图7 滑膜面和signum模块的时域响应 图8 1s左右时的滑模面和signum响应的局部放大图如图9所示,系统大约经过12个仿真步到达1s处时,步长急剧缩小至接近于零。这样,由于仿真步长太小,系统便在不连续处形成了过多的点,超出了系统可用的内存和资源,使得系统进展缓慢,仿真停滞不前。 图9 系统仿真步长3 问题的对策通过以上分析,我们总结

20、、概括得出了以下系统仿真速度慢的原因:(1 系统方程中存在不连续函数sgn(s;(2 Simulink仿真模型中存在能够产生过零通知的signum模块;(3 采用的变步长求解器具有过零检测并自动调整步长的功能。3.1 解决策略因此,基于以上原因,我们提出了以下4种解决策略:(1 取消Zero crossing detection功能;(2 采用fixed-step求解器;(3 采用不能够产生过零通知的Fcn函数模块;(4 柔化sgn(s函数,使其连续化。以上所提出的4种解决策略中,采用单独任何一种或几种都可行。以下,我们尝试采用第3种策略来观察系统的仿真效果。即采用Fcn函数模块来代替图3中的

21、其他模块所表示的数学关系,搭建的仿真模型如图10所示: 图10 采用Fcn函数模块后的系统仿真模型对于函数fun,我们取2x,1x,a分别作为输入u1, u2, u3。根据式(2可得函数fun的表达式为:-9.8*u2+u3对于函数fun1,我们取2x,1x,1c,分别作为输入u1, u2, u3, u4。根据式(7可得函数fun1的表达式为: 9.8*u2-u3*u1-u4*sgn(u3*u2+u1取11c=,1=,系统仿真结果如图11所示: 图11 系统时域图2006年7月 系 统 仿 真 学 报 July, 2006由于Fcn 模块不支持过零,所以,系统在不连续的情况下仍然能迅速完成仿真

22、。除此之外,我们还采用了fixed-step 求解器和变步长下置Zero crossing detection 为off 的仿真方法,其仿真结果亦如图11所示。可见,这前三种方法的仿真结果一样。由上面的仿真结果可以看出:(1 取消系统的过零检测功能之后,仿真速度快; (2 但是由于仍然存在不连续模块signum ,所以加速度存在较大的抖振。3.2 sgn(s函数连续化 由于加速度存在剧烈的抖振,使得控制系统难于工程实现,为了能消除系统因为不连续性而存在的抖振现象,同时加快仿真速度,且不影响系统的仿真效果,我们设计了第四种控制策略-sgn(s函数连续化。连续化前后的函数sgn(s函数和sat(s

23、函数如图12所示: 图12 sgn(s和sat(s函数其中,sat(s函数的表达式为:1(1s sat s s = <(8 为大于零的正数,且。当取无穷小时,sat(s 函数便非常逼近sgn(s函数。采用sat(s函数来代替sgn(s函数后的系统仿真模型如图13所示: 图13 柔化sgn(s函数后的系统仿真模型根据式(8可得函数fun 的表达式为:sin(u/0.1414。其表示,s 时sat(s的值。仿真结果为:(I 当0=时,连续函数sat(s即转化为符号函数sgn(s,仿真时将出现如前面所述的问题(仿真速度慢、抖振等。(II 当0.1=时,其它参数与图3中相同,可得系统仿真结果如图

24、14所示:图14 采用连续函数(0.1=后的系统时域图仿真结果说明:采用连续函数sat(s代替不连续函数sgn(s后,系统仿真速度加快,不影响系统控制效果,且加速度不存在抖振现象,易于工程实现,满足我们所期望的结果。(III 当0.0001=时,系统的仿真结果如图15所示:图15 0.0001=时的系统仿真结果仿真结果说明:当0.0001=时,加速度出现剧烈的抖振,性能接近理想的滑模变结构控制。通过以上的仿真结果,我们可以总结如下:(1 采用sgn(s函数连续化的方法,系统仿真速度快; (2 越小,系统仿真结果越逼近原理想的滑模变结构控制系统,但抖振也越大,相应的仿真速度亦越慢;反之,越大,仿

25、真结果失真也越大,但系统抖振越小,相应的仿真速度亦越快。因此,我们建议:(1 如果希望既不影响系统的仿真速度,又不损失精度,0.0001为宜;第 18 卷第 7 期 2006 年 7 月 Vol. 18 No. 7 系 统 仿 真 学 报 表4 步长 h=0.1 求解例 2 在 x July, 2006 4 数值试验 例 1（参见5） 常数延迟微分方程初值问题 = 10 处的整体误差. x = 10 处整体误差 1.6235864E-4 6.70327E-4 y '(x = 1000y(x + 997e y(x 1 +1000 997e 3x y(x =1+ e 3 3 0 x4 x0

26、 h=0.2 上面两个例子验证了我们所构造的方法的收敛阶， 并且 验证了方法在求解刚性问题时是有效的。 问题的精确解是 y ( x = 1 + e 的整体误差列在表 3 中。 表3 步长 h=0.1 h=0.01 3 x 。 用表 2 列出三阶方法求解例 1，在积分区间0,4内得到 参考文献: 1 H Podhaisky, B A Schmitt, R Weiner. Design, analysis and testing of some parallel two-step W-methods for stiff system J. Applied Numerical Mathematics

27、(S0168-9274, 2002, 42: 381-395. 2 3 H. Podhaisky, R. Weiner, B.A. Schmitt, Two-step W-methods for stiff ODE systems J. Vietnam J. Math. (S0866-7179, 2002, 30: 591-603. Bernhard A. Schmitt , Rüdiger Weiner, Parallel Two-Step W-Methods with 4 Peer Variables J. SIAM Journal on Numerical Analysis(S

28、0036-1429, 2004, 42: 265-282. K J int Hout, M N Spijker. Stability analysis of numerical methods for delay differential equations J. Numer. Math. (S0006-3835, 1991, 59: 807-814. 5 6 曹学年，刘德贵，李寿佛. 求解延迟微分方程的 ROSENBROCK 方法的渐近稳定性J. 2002，14(3：290-292. E Hairer, S P NØrsett, G.Wanner. Solving Ordinary

29、 Differential Equations I. Nonstiff problems M. Berlin: Springer-Verlag , 1993. 求解例 1 在区间0,4上的整体误差 区间0,4上的整体误差 1.1762236E-05 1.5437829E-08 例 2（参见6） y1' ( x = y1 ( x y2 ( x 1 + y2 ( x 10 ' y2 ( x = y1 ( x y2 ( x 1 y2 ( x ' y3 ( x = y2 ( x y2 ( x 10 100 y3 ( x 当 x 0 时， y1(x = 5 ， y2 ( x = 0.1 ， y3 ( x = 1 ，积分区间 是 0,10 ，我们用表 1 列出 GP-稳定的二阶方法求解例 2。 例 2 是文献6的 348 页中传染病模型的例子， 我们在这 里增加了-100y3(x这一项，使问题成为刚性的。表 4 列出了 例 2 的计算结果。 （上接第 1754 页） (