第8章 控制系统仿真工具

第8章控制系统仿真工具8.1 仿真工具概述计算机是系统仿真最重要的工具，用于仿真的计算机基本可分为通用计算机和专用计算机。

常用的仿真语言有C语言、Visual Basic和Matlab。

C语言主要用于仿真研究后直接进行实际工程应用(单片机、工控机等)且对界面要求不高的情况；Visual Basic主要用于仿真研究后直接进行实际工程应用(微机、工控机)且要求界面设计美观的情况；而Matlab主要用于纯仿真研究，以便得出一些具有参考意义的结论。

除计算机以外，仿真工具还包括两类仿真器：专用的物理仿真器和用于训练操纵人员的所谓操纵训练仿真器。

8.2 VB语言仿真(1) Visual Basic的优点1976年, 现任美国微软公司总裁Bill Gates(比尔·盖茨)先生和Paul Allen开发出了Basic 语言的早期版本。

1998年12月微软公司正式推出Visual Basic 6.0中文版本，尤其是近3年，Visual Basic几乎是两年更新一个版本，说明Visual Basic程序开发的投资力度和受重视的程度。

可以说Visual Basic是Bill Gates先生的招牌，他是靠Basic语言起家的，故其发展的速度和投资力度是有保证的；另一方面，Visual Basic的迅猛发展也是其应用广泛和受用户喜爱的一个标志。

用VB编程具有以下5个优点：1)用户图形界面编程工作量少。

一个实际的软件编程可分为图形界面编程和非图形界面编程两部分，分别代表用户可见的和不可见的部分。

由于VB编程面向用户图形界面，用户只需按设计风格和实际要求将VB系统环境所提供的基本元素(包括按钮、图片、文本框、标签等)放在窗体中，然后更改其属性表即可，就像用各种事先做好的模具来绘画一样，只需按脑海中构想的图找到相应的模板在图画纸的相应位置上画好，然后涂上颜色就行了。

所以说用VB编程，用户可见部分的图形界面编程的工作量大大降低，从而缩短了软件的开发周期。

2)VB编程以事件驱动为机制。

用户的任何一个动作都产生一个事件，如按鼠标左键、击键等，如果给一个事件编程了，则当该事件发生时就会产生执行效果。

可以给所有可能发生的事件都编程，则任何一个事件发生都会产生效果。

就像玩捅开纸格中奖游戏一样，如果每一个纸格后面都藏有奖品，则每捅开一次纸格都会中奖，只不过奖品不同而已。

也有可能纸格后面没有藏奖品(相当于某一事件没有编程)，则捅开该纸格后将得不到奖品(相当于该事件发生后将没有执行效果)。

同时其事件驱动的编程机制使得程序设计起来容易，读起来直观，使用起来方便。

3)VB采用面向对象的程序设计方法，程序代码具有良好的可重用性、可扩充性和可修改性。

4)VB 6.0作为Windows编程技术之一，具有Windows环境所具有的5大优点：即标准的图形用户界面、动态链接(DDL)、多任务、设备独立性及直接操作特性。

5)由于VB的用户越来越多，很多计算机软件公司相继开发出适于VB的工具箱，目前由第三方(即非Microsoft公司和用户以外的第三计算机软件开发方)开发的VB工具箱已达上百种，故VB程序设计随时间的推移会越来越简单。

关于Visual Basic 程序设计的学习请参见文献[28]和文献[29]。

(2)用VB 进行计算机仿真的工作步骤1)将要仿真的系统数学模型写成微分方程组形式，并编写成函数或普通过程。

2)选择适当的系统离散化方法，并编写成函数或普通过程。

3)计算或估算出系统的采样步长大小。

4)对控制器选择适当的离散化方法，写出递推公式。

5)给定输入，按反馈系统信号传递方式，由递推公式，求出输出。

6)判断输出的精度是否满足要求，如果满足要求转步骤7，否则转步骤5。

7)画仿真曲线，输出仿真数据，并根据仿真结果进行分析，给出结论。

在进行连续系统数字仿真时，离散方法的选择十分重要。

但是究竟如何选择，至今没有一种具体办法。

所以，只能根据实际系统的性能、所要求的精度和计算机的速度、内存等性能等实际情况的不同具体选择一个较为合适的离散方法。

下面给出4阶龙格-库塔法的VB 程序清单：Sub rk4 ()Dim i, j As Integeru(0) = .5 * hu(1) = u(0)u(2) = hu(3) = hu(4) = u(0)For i = 0 To n '赋数组初始值f(i) = 0#xr(i) = x(i)yw(i) = x(i)Next iFor i = 0 To 3fss '计算k1,k2,k3,k4，fss 为数学模型微分方程组过程名For j = 0 To nxr(j) = yw(j) + u(i) * f(j) '为计算k2,k3,k4作准备x(j) = x(j) + u(i + 1) * f(j) / 3# 'x(j)为龙格-库塔解Next jNext iEnd Sub说明：该4阶龙格-库塔VB 普通过程代码适用于对具有由n 个微分方程组成的数学模型进行离散求解，是一个通用的程序，读者一步一步验算可知其公式与前面给出的4阶龙格-库塔公式是完全一致的。

虽然读起来不那么直观，但程序中应用了一些编程技巧和经验，所用语句、变量都最少，具有节省内存、运行速度快的优点。

当然要想正确应用该程序，还需：≠ 将h, n, u(4)定义成模块级或全局的常数和数组，将f(微分方程的微分变量数组), xr(微分方程的状态变量数组), yw(中间变量数组), x(初值变量及解变量数组)定义成具有n 维的模块级或全局数组，注意n 必须为具体的数。

≡ 在名为fss 的普通过程中按标准状态空间方程的形式即BU AX X+=&编写数学模型的微分方程，f 数组代表微分变量数组，xr 数组代表状态变量数组。

≈ 在事件过程中编写用龙格-库塔方法解算数学模型的计算机仿真程序，包括赋变量初值、循环解算并判断解是否满足精度要求、解的数据及图形曲线输出等。

【例8.2.1】 用VB 6.0对船舶航向自动舵控制系统进行计算机仿真，目的是测试所设计的航向自动舵控制算法是否满足所要求的性能指标，假设其控制算法采用PI 控制，而船舶运动仿真数学模型采用第6章的三自由度模型。

控制精度要求：航向误差小于0.1︒，上升时间小于250 s ，最大超调量不大于40%。

系统框图如图8.2.1所示。

仿真要求给出航向fai 及舵令delt 的仿真曲线结果，并分析仿真结果。

图8.2.1 船舶航向自动舵控制系统框图控制算法采用：PI 控制，即sk k K i p +=，其中比例系数和积分系数根据控制的仿真效果自己调整，控制器的离散建议用Tustin 变换，利用其简单、稳定的优点。

船舶运动仿真数学模型采用三自由度模型，离散算法采用RK4或RK5，利用其适用于非线性模型的优点。

步长自己按经验公式计算。

数学模型的数据来源于大连远洋运输公司的“阳澄湖”号油船的数据，全部采用状态空间数学模型，为了简化模型，暂不考虑风、浪、流等干扰项，即使用名义数学模型，现给出其模型参数如下：)3()3()3()3()3()3(.X C B X A X =+=fai delt (8-2-1) 其中T )3(][fai rv =X r v ,分别为船舶的横漂速度和艏摇角速率⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡----=0100089.0001.00846.1027.0)3(A ⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-=0002.0042.0)3(B []100)3(=C 用VB 6.0编程实现的步骤如下：1)计算采样时间步长，根据第7章给出的经验公式，首先绘制出船舶运动数学模型的Bode 图如图8.2.2所示，求出03.0=c ωrad/s ，则采样时间步长h 取6 s 。

2)采用Tustin 变换将K 离散化得fai10-4 10 -3 10 -2 10 -1 10 0 10 1Frequency / (rad ⋅s -1) G a i n /d B 10 -4 10 -3 10 -2 10 -1 100 10 1Frequency / ( rad ⋅s -1) P h a s e /(︒)图8.2.2 船舶运动数学模型的Bode 图⎪⎪⎭⎫ ⎝⎛+-=--+==11112z z h s ip s k k dfai delt K)1()2/()()2/()1()(---++-=k dfai h k k k dfai h k k k delt k delt i p i p3)确定仿真结束条件，本例通过判断仿真时间长短作为仿真结束条件，要求仿真2 400 s 即400步，系统仿真达到稳态后为直线，观察价值不太，如果在2 400 s 内系统还未达到稳态，则认为控制器参数不好，应重新调整PI 控制器的比例系数和积分系数。

4)启动VB 6.0，更改窗体的Caption 属性为“船舶航向自动舵控制系统仿真”，Borderstyle 属性改为1-Fixed Single ，并将窗体拖至足够大用以显示两个曲线图。

5)打开代码窗口，编写代码如下：'在general 中定义模块级变量Const h! = 6# '采样时间为6.0 sDim u(4) As SingleConst n% = 2 '3个微分方程Dim f(2), xr(2), yw(2), x(3) As DoubleDim delt, dfai, fair, fai As SingleDim aa, bb As IntegerSub fss() '数学模型普通过程f(0) = -0.027 * xr(0) - 1.846 * xr(1) + 0.042 * deltf(1) = -0.001 * xr(0) - 0.089 * xr(1) - 0.002 * deltf(2) = xr(1)End SubSub rk4() '4阶龙格-库塔普通过程Dim i, j As Integeru(0) = 0.5 * hu(1) = u(0)u(2) = hu(3) = hu(4) = u(0)For i = 0 To n '赋数组初始值f(i) = 0#xr(i) = x(i)yw(i) = x(i)Next iFor i = 0 To 3fss '计算k1,k2,k3,k4，fss为数学模型微分方程组过程名For j = 0 To nxr(j) = yw(j) + u(i) * f(j) '为计算k2,k3,k4作准备x(j) = x(j) + u(i + 1) * f(j) / 3# 'x(j)为龙格 库塔解Next jNext iEnd Sub'绘图子程序Sub ht0(zsx As Integer, zsy As Integer, yxx As Integer, yxy As Integer, xi As Integer, yi As Integer, x0 As Single, xstep As Single, page As Integer, y0 As Single, ystep As Single, xtx As Integer, xt As String, yty As Integer, yt As String, titlex As Integer, title As String)'zsx左上角x坐标,zsy左上角y坐标,yxx右下角x坐标,yxy右下角y坐标,xi为x轴标度坐标间隔,yi为y轴标度坐标间隔,x0为x轴标度数据起点，xstep为x轴标度数据步长,page为页号(曲线连续多页显示时用)，y0为y轴标度数据起点，ystep为y轴标度数据步长,xtx为x轴说明起始坐标,xt为x轴说明字符串，yty为y轴说明起始坐标,yt为y轴说明字符串,titlex为图的标题起始x坐标,title为图的标题Dim i, j, hxw, hxt, sxw, sxt, titley, xty, ytx, pagex, pagey, dataxx, dataxy, datayx, datayy As Integerhxw = 20 '横虚线宽hxt = 40 '横虚线加上虚线缝隙总宽sxw = 1 '竖虚线宽sxt = 2 '竖虚线加上缝隙总宽titley = 8 '标题y坐标距上边框距离xty = 2 'x轴说明坐标距底边框距离pagex = yxx '页码x坐标pagey = zsy '页码y坐标ytx = 200 'y轴说明坐标距左边框距离dataxx = 100 'x轴数据标注间隔dataxy = 2 'x轴数据标注距底边框距离datayx = 120 'y轴数据标注距左边框距离datayy = 2 'y轴数据标注间隔(zsx, zsy)-(yxx, yxy), RGB(0, 255, 0), B '画曲线边框For i = 0 To ((yxy - zsy) / yi - 2) '画横的虚线网格For j = zsx To yxx - hxw Step hxt(j, zsy + (i + 1) * yi)-(j + hxw, zsy + (i + 1) * yi), RGB(0, 255, 0)Next jNext iFor j = 0 To ((yxx - zsx) / xi - 2) '画竖的虚线网格For i = zsy To yxy - sxw Step sxt(zsx + xi * (j + 1), i)-(zsx + xi * (j + 1), i + sxw), RGB(0, 255, 0)Next iNext jntX = titlex '标题显示ntY = zsy + titleytitlentX = xtx 'x轴说明ntY = yxy + xtyxti = 0Do While i <= (yxy - zsy) / yi '标注y轴数据ntX = zsx - datayxntY = yxy - i * yi - datayyStr$(y0 + i * ystep)i = i + 1Loopj = 0Do While j <= (yxx - zsx) / xi '标注x轴数据ntX = zsx + j * xi - dataxxntY = yxy + dataxyStr$(x0 + j * xstep)j = j + 1LoopntX = zsx - ytx 'y轴说明ntY = ytyytEnd SubSub ht(a As Double, b As Double, i As Integer) '动态绘图If i = 1 Then '准备航向误差初始坐标值CurrentX = 300CurrentY = 30Else '准备航向误差当前坐标值CurrentX = (i - 1) * h + 300CurrentY = -aa + 30End IfLine -(i * h + 300, -CInt(a) + 30), RGB(255, 0, 255)If i = 1 Then '准备舵令初始坐标值CurrentX = 300CurrentY = 90Else '准备舵令当前坐标值CurrentX = (i - 1) * h + 300CurrentY = -bb + 90End IfLine -(i * h + 300, -CInt(b) + 90), RGB(0, 0, 255)aa = CInt(a) '当前坐标更新bb = CInt(b)End SubPrivate Sub Form_Click()Dim delt1, dfai1, kp, ki As SingleDim i As Integerkp = -0.1 'PI控制器参数设置，负号是因假设舵令方向与实际方向相反ki = -0.0001olor = RGB(255, 255, 255) '设背景色为白色edraw = True '图形自动重画Scale (0, 0)-(3000, 120) '定义窗体坐标系'绘制航向误差曲线ht0 300, 10, 2700, 50, 200, 10, 0, 200, 1, -20, 10, 2800, "(秒)", 2, "dfai(度)", 1300, "航向误差曲线"'绘制舵令曲线ht0 300, 70, 2700, 110, 200, 10, 0, 200, 1, -20, 10, 2800, "(秒)", 62, "delt(度)", 1500, "舵令曲线"fair = 10 '设定航向假设为10度fai = 0 '设初始航向为0度delt = 0 '设初始舵令为0度delt1 = 0dfai1 = 0'控制器输出递推解算For i = 1 To 400 '仿真400步即2 400 sdfai = fair - faidelt = delt1 + (kp + ki * h / 2#) * dfai - (kp - ki * h / 2#) * dfai1delt1 = delt '递推变量更新dfai1 = dfairk4 '解算船舶运动数学模型fai = x(2)x(3) = deltht x(2), x(3), i '绘制航向误差及舵令曲线Next iEnd Sub6)按F5运行程序，点一下窗体，仿真曲线如图8.2.3所示。