第四章 面向结构图的数字仿真法

1、第四章 面向结构图的数字仿真法对一个控制系统进行研究,其中一个很重要的问题就是考察系统中一些参数 改变对系统动态性能的影响,面向微分方程的仿真方法很难得到这一点。这主要是由小回路的传递函数 得到的全系统大回路的传递函数之间的参数对应关系将变得非常复杂。其次,将复杂系统中诸多小回路 化简求出总的系统模型也是十分麻烦的,更何况对于非线性系统,或难以用非数学模型描述的系统,则 无法找到系统的总的闭环模型。本章介绍两种由一些典型环节构成的复杂系统仿真的方法。在这类仿真程序 中,先将仿真这些典型环节特性的仿真子程序编制好;在仿真时,只要输入各典型环节的参数以及环节 间的连接关系的参数便可以作系统的仿真。

2、这就是面向结构图的数字仿真法,它可以解决上述困难,且 具有一些优点:(1很容易改变某些参数环节,便于研究各环节参数对 系统的影响。(2不需要计算出总的传递函数,并且可以直接得到各 个环节的动态性能。(3系统中含有非线性环节时也比较容易处理。本章第一节介绍面向结构图仿真各典型环节仿真模型的确定。第二节介绍面 向结构图模型离散相似法仿真的方法。第三节介绍对于含有典型非线性环节的处理方法。第四节介绍连 续系统结构图仿真方法、程序的编制及应用。4.1 典型环节仿真模型的确定在第二章第四节中已经介绍了状态方程离散化的方法,即对一个状态方程加入虚拟的采样器和保持器,当采样频率合适时则可实现信号重构。面向结

3、构图仿真方法其基本思想就是将结构图化简为各个典型环节组成,然后在各个典型环节前加入虚拟的采样器和保持器使各环节独自构成一个便于计算机仿真的差分方程。本节就是求出个典型环节对应 的离散状态方程得系数矩阵,即。1.积分环节 其中。离散状态方程为 2.比例积分环节 其中: 3.惯性环节 其中:环节的状态方程和输出方程为 离散状态方程为 4.比例惯性环节 其中:。状态方程为 显见状态方程与惯性环节一样,故的计算也一样,仅输出方程不一样,故得离散状态方程为 可见高阶环节均可用前述几种典型环节获得。4.2 结构图离散相似法仿真面向结构图模型的离散相似法仿真除了需要建立典型环节的差分式外,还需要建立能描述系

4、统连接方式的方程。在上一节的基础上,本节将进一步介绍系统连接矩阵的建立和面向结构图模型的离散相似法仿真方法以及计算程序的实现。一、连接矩阵上一节,介绍了环节离散化方法以及所得到的差分方程模型的形式。但这仅仅表示了各个单独环节输入和输出之间的关系。为了实现面向结构图离散相似法仿真,还必须把这些环节按照系统结构图的要求连接起来、以保证正确的计算次序。设系统的第i各环节输入、输出分别用表示,y0为系统的外部输入量,则式中W是一个n×(n+1维长方矩阵。这是把表示输入信号与系统连接情况的W0矩阵放在原连接矩阵的第一列,也就是 W ij表示第j个环节输入之间的连接方式。所示。如果已知各环节的传

5、递函数,则很容易将其离散化,各个环节的输入输出关系为 二、仿真程序的设计把不同类型环节的离散系数的计算分别编成子程序。在程序中引入一个标志参数,表示该典型环节的类型,假设一个通用程序只包括下列四种典型环节,且与典型环节对应关系如下当=0时,表示第I个环节为积分环节。当=1时,表示第I个环节为比例积分环节。当=2时,表示第I个环节为惯性环节。当=3时,表示第I个环节为比例惯性环节。由前述可知,对于=0和=1两种典型环节,计算状态变量x的公式相同,只 是它们输出变量计算公式不同。而同样对于=2和=3的典型环节,也是计算状态变量x的公式相同,仅仅是输出方程不同。在步长取定后,典型环节的离散系数,就仅

6、是典型环节的参数(时间常数、放大增益的函数,可以预先根据典型环节的类型分别编成子程序,仿真时即可根据方便地调用。系统的连接情况,仍用连接矩阵W来描述。 按系统的典型环节离散化仿真,其主要优点是:(1各个环节的离散状态方程系数计算简单,而且可以一步求出,不像龙格-库塔法那样,每一步都要重新计算龙格-库塔系数,因而计算量相对来说较小。(2由于各个环节的输入量u i,输出量y i每一步都可求出,所以很容易推广到含有非线性环节的系统仿真中去。该方法的主要缺点是计算精度低。因为每个环节的输入实际上都是使用了它们的近似值(举行近似或梯形近似,故仅有一阶或二阶精度,这会带来计算误差,而且环节越多,误差越大。

7、这一点下面还将进一步分析。另外,需要指出的是,当输入采用梯形近似法时,需要用到u(n+1来求取u(n,u(n=u(n+1-u(n/T,这通常是难以办到的。于是在仿真中有时只得采用简单的向后 差分的方法来计算,即。由于本来的定义是表示在n T(n+1T 区间输入信号的平均变化速度,所以用向后差分的方法来计算,实际上使用前一个周期(n-1Tn T 的输入信号的平均变化速度来近似代替周期n T(n+1T的输入信号变化速度,相差一个采样周期。这显然会使计算误差增大。三、仿真算例及分析 根据经验公式,大约可达到0.5%左右的精度,为系统开环频率特性的剪切频率。在此例中因此可选T=0.025s。输入数据有

8、 连接矩阵 仿真参数采样周期仿真时间打印、显示时间间隔输出环节号10.01 10 1 1输入以上三组参数后,便可在计算机上仿真。四、采用补偿器提高模型精度和稳定性的方法系统的离散化过程,就是在连续系统中加入虚拟的采样开关和保持器。由于保持器不可能完整无误地将连续信号重构出来,因此必然会产生仿真误差。一般来讲,采样间隔越大,仿真的误差也就越大。为了减少误差,很自然地就想到是否能在这个仿真器中模型中加进校正补偿环节。一般所加入的补偿器应尽可能好地抵消经过采样-保持器所造成的失 下面以积分环节为例来说明这种方法的基本原理。 对做一次近似,即去 写成差分方程 选择不同的,可得各种不同数值及分公式。比如

9、 所示,则有 这就是亚当斯公式。 当将这种方法用于复杂系统时,为获得仿真模型,其基本步骤是: (2 求出该图所示的离散化系统的脉冲传递函数,并列出它的差分方程,这就是仿真模型。(3 用高阶的龙格-库塔法计算该系统的相应,将它作为一个标准解,然后给出不同的 ,计算仿真模型的响应,并将它与标准解进行比较,直到误差达到最小为止。利用上述步骤仅仅是计算出了系统的输入量y。如果不仅对y感兴趣,而且对于系统中的其他变量也有兴趣,那么就必须将系统分成几个部分,每部分都要加虚拟的采样器及保持器。至于校正补偿环节则按一般系统的校正原则,可以对每一个小闭环加一个。调整时,一般是先调外环的。调整的目标是要求所获得的

10、仿真模型在较大的计算步距时仍能 最好地与实际模型相接近。地选取可采用第六章所介绍的参数寻优程序来确定。4.3 非线性系统的数字仿真在本章的上一节中曾提到,利用离散相似法编制的仿真程序虽然精度低,但是却可以十分方便地推广应用到这类非线性系统中去,其主要原因是在仿真计算程序中,每走一步,各个环节的输入量及输出量都将重新计算一次。因此非线性环节子程序很容易加入到仿真程序中去。下面首先介绍典型的非线性环节的仿真。一、非线性环节仿真子程序1.饱和非线性 2.失灵区非线性 3.齿轮间隙(磁滞回环非线性(图 设为上一次的输入,为上一次的输出。若当 >0时,且 C1,则 = C1 即,若只满足前一个条件

11、,而不满足后一个条件,则是工作在由左边的特性过渡到右边的特性上。 若当 <0时,且+C1,则 = + C1 二、含有非线性环节的离散相似法仿真程序的计算方法当系统中有上述典型环节时,本章第二节讲的离散相似法仿真程序要作如下修改:(1 对每个环节要增设一个参数,表示第个环节的入口或出口有那种类型的非线性环节。(2 对每个环节要增设一个参数,表示第个环节的入口的那个非线性环节的参数C i,当第个环节入口没有非线性时,=0。因此在输入数据时,对于每一个非线性环节都要同时送 8个数据。其中为线性环节 的系数,的含义为: =0 表示该环节前、后无非线性环节。 =1 表示该环节前有饱和非线性环节。

12、=2 表示该环节前有失灵区非线性环节。 =3 表示该环节前有齿轮间隙非线性环节。 =4 表示该环节后有饱和非线性环节。=5 表示该环节后有失灵区非线性环节。=6 表示该环节后有齿轮间隙非线性环节。三、非线性系统仿真举例第一步:确定系统各个环节号。本例除第一个环节前有非线性环节外,其余都为线性环节。 第三步:运行程序根据提示输入数据。输入的数据有:(1各环节参数,按第二种情况考虑,非线性环节为饱和特性,即 u0(i y0(i Z(iD(iC(iA(i B(i0.1 1 0.5 1 0 0 10 1 1 0 0 0 02 1 2 0 0 0 010 1 10 0 0 0 0(2输入连接矩阵数据,即

13、 (3输入仿真参数。根据经验公式,采样周期按各环节最小的时间常数的1/10选取,本例中最小时间常数为0.1,故采样周期选为0.01。仿真时间去10s,且观察第1号、4号环节输出,因此仿真参数输入时如下:采样周期仿真时间打印、显示时间间隔输出环节号0.01 10 1 1,4第五步:结果分析。将以上三组数据输入到仿真程序中,运行后可得到数据结果。通过该程序的仿真结果,可以分析系统中典型非线性环节对系统的影响。1.饱和非线性对系统过渡过程的影响当自动控制系统中(非条件稳定系统存在饱和元件时,此时系统的稳定性将变好,而快速性将变坏,也即超调量将减小,而过渡过程时间增加。这与自动控制原理理论分析结果相同

14、。2.失灵区非线性对系统过渡过程的影响根据调解原理分析可知:若系统中具有失灵区非线性环节,那么系统的动态品质将变坏,而对稳定性影响不大。其原因是:(1 由于有失灵区,在过渡过程的起始段,相当于减小了系统的开环增益,故过度过程变缓。(2 当输入量接近稳定时,放大器处于稳定区,系统处于失控状态,控制作用为零,故超调量将略微增大。(3 由于放大器有失灵区,故在过渡过程中有尾部,系统也处于失控状态,因此将出现一个很长的尾巴。即从系统进入失灵区到输出量进入稳态值区(±5%,输出量变化十分缓慢。3.齿轮间隙非线性对系统的影响由于存在齿轮间隙,当系统的输出值超过稳态值时,因系统有反向调节的趋势,输

15、出将维持不变一直要等非线性环节的输入走完间隙时输出才能下降。而当输出值反向偏离稳态值时,系统邮政向调节的趋势,输出又将维持不变,一直要等非线性环节的输入走完间隙时输出才能回升。其结果,系统将会在稳态值附近以某一幅度和频率进行震荡,即系统始终在一个极限环内运动,而无法稳定下来。4.4 连续系统的结构图仿真及程序本节将介绍一个面向结构图的数字仿真程序包CSS(Continuons System Simulation。该程序包是1981年5月引入我国的。该程序包首先被移植在国产的DJS-130机上,后因BASIC 语言运行速度太慢,国内有关单位将该程序包翻译成FORTRAN语言,并对原程序作了大量的

16、修改和补充,取名为CSSF,ZFX等。CSSF仿真程序包的主要特点:(1 配备了多种积分方法。例如,定步长龙格-库塔法、边步长梅森法、定步长汉明法等。(2 可以面向多种形式的数学模型。例如,状态方程,一阶及二阶传递函数,n阶传递函数或微分方程,对所有这些原始方程,只需直接输入系统,不必进行交换。(3 增加了延迟、微分、一元及二元函数发生器等。(4 增加了采样系统仿真功能,并且有参数优化模块。(5 具有非线性两点边值自动求解功能。CSSF程序包提供了约40种不同类型的运算块,包括模拟机中的积分器、比例器、三角函数、对数、乘法、除法以及各种非线性函数运算模块,故一套CSSF程序包的计算功能相当于一

17、台的大型模拟机。除此之外CSSF程序包还包括几块用户自定义块意识和用户的特殊需要。用CSSF进行仿真的最大优点就是使用方便。用户只需将仿真系统转化成程序包所含有的运算功能块所组成的仿真结构图,启动程序包,按指定的方式输入个功能块的编号、类型参数、连接方式,以及积分步长、仿真时间等参数,就可以得到仿真结果。可以说,CSSF仿真程序是一种功能全、性能高的程序包,可用于大型、复杂的连续或采样系统的仿真,但该程序包的结构复杂。Micro-CSS(MCSS是CSS的微型化,它是用BASIC语言编写的。基本的MCSS约有80条语句,而扩展MCSS约有200条语句,并采用C语言编写。它包括了CSS程序包的主

18、要内容,并且参数的输入风格也与CSS相同,故通过MCSS程序的学习,也可解剖CSS的结构,并掌握CSS的使用方法。 下面结合具体例子说明构成MCSS仿真图的方法、系统参数、连接情况,以及仿真参数的输入方式。 初值:。试画出仿真结构图,并写出仿真数据。解 1.仿真图构成方法由于原微分方程可改写为 2.各功能块连接方式输入数据各功能块连接方式输入数据的格式为DATA I,16,J,K,L其中I-功能号顺序号;16-功能号类型号;J-同第1输入端(C2连接的模块号;K-同第2输入端(C3连接的模块号;L-同第3输入端(C4连接的模块号。这部分数据输入结束时,用全零行表示,即DATA 0,0,0,0,

19、0本例的数据格式为10 DATA 1,1,0,0,012 DATA 2,3,1,3,214 DATA 3,3,2,0,016 DATA 0,0,0,0,0 3.参数部分数据输入参数部分数据输入的格式为 其中:I功能块顺序号 这部分数据输入结束时,同样用全零行表示。注意:积分块(3型功能块的第一参数输入为初始条件数据。当某一功能块中p1,P2, P3三个参量都没有时,可以不输入该块的数据。例如,正弦函数块、余弦函数块、乘法器等。本例的数据格式为60 DATA 1,1,0,062 DATA 2,0,0.1, 164 DATA 3,0,0,066 DATA 0,0,0,04. 仿真控制数据输入仿真控

20、制数据输入的格式如下:80 DATA D,T1,T2,K1,K2,K3,K4其中D积分步长;T1仿真时间;T2打印、显示时间间隔;K1K4输出量所在模块。 本例中,如果设计步长为0.1s,仿真时间为10s,每秒显示一次,输出量为,则数据格式为90 DATA 0.1,1,3,2,0,0其输出结果如下:TIME OUTPUT3 OUTPUT2 OUTPUT0 OUTPUT00 0 0 0 00.999 996 0.365 040 3 0.721 808 0 01.999 999 1.103 401 0.811 573 8 0 02.999 998 1.932 805 0.830 502 8 0 0

21、3.999 997 2.743 704 0.785 397 5 0 04.999 996 3.496 954 0.719 563 2 0 05.999 995 4.181 771 0.650 230 8 0 06.999 994 4.798 479 0.583 958 2 0 07.999 993 5.351 465 0.522 966 2 0 09.000 016 5.846 339 0.467 740 8 0 0 数据文件为10 DATA 1,1,0,0,012 DATA 2,2,1,8,114 DATA 3,16,2,0,016 DATA 4,3,0,3,018 DATA 5,3,4,5

22、,020 DATA 6,3,5,6,022 DATA 7,3,0,6,724 DATA 8,2,6,7,026 DATA 0,0,0,0,050 DATA 1,1,0,052 DATA 2,1,-1,054 DATA 3,1,-1,056 DATA 4,0,4,058 DATA 5,0,-2,060 DATA 6,0,-10,062 DATA 7,0,0.4,-0.164 DATA 8,1,1,066 DATA 0,0,0,0,090 DATA 0.01,5,0.5,8,0,0,0其输出结果为TIME OUTPUT3 OUTPUT2 OUTPUT0 OUTPUT00 0 0 0 00.499

23、999 8 2.753 746E-02 0 0 00.999 999 3 0.110 005 2 0 0 01.499 999 0.226 193 0 0 01.999 998 0.366 024 2 0 0 02.499 998 0.525 232 6 0 0 02.999 998 0.701 728 2 0 0 03.499 997 1.101 851 0 0 04.500 008 1.323 749 0 0 0为了方便读者进一步该程序的使用方法和应用范围,并熟悉变成技巧,本节对该程序的组成及功能作更深入的分析和说明。MCSS程序共有以下四部分组成:(1 输入数据块;(2 自动安装模块顺序

24、块;(3 运行程序快;(4 输出打印块;下面分别说明以上四部分程序。1.输入数据块这部分主要有读语句组成,输入的数据有三部分,按先后顺序为(1 连接方式数据:I,16,J,K,L(2 典型块参数数据:I,y(1,y(2,y(3(3 仿真控制数据:D,T1,T2,K1,K2,K3,K4其中(1,(2两类输入数据结束标志为零。数据段可放在程序的前部或后部,但三部分数据的次序不能改变。类型数目编号延时块N4ID(N4积分块N9V(N9常数块N8N(N82.自动安排各模块计算顺序这一段程序的主要功能是编排各模块计算顺序并赋予顺序号。用户在完成仿真结构图后,先要给你一个计算顺序,即编号。至于这一顺序号是

25、否合理,程序会自动进行判断,并作出合理安排。前排的原则是,该模块的三个输入信号是否已有确定值,若无,则这一模块必须排在以后计算。但有三种类型块例外:常数块(1、积分块(3、和延迟器(4。它们的输出值可由给定的参数和初始值来确定,不必依赖于输入值。 在进入这段程序以前,除了常值信号源外(因为它们不必计算,故没有安排计算顺序的必要,其他运算块的类型代码(程序中用U(I,1表示,均为负值,以表示尚未安排好的顺 下面用一个例子来说明MCSS程序中各模块计算顺序自动安排的过程。 设初始条件为 试分析计算机排序的方法和结果。解:首先将方程改写为 程序先按用户编号1号开始检查,因为1号比例器只有第3输入端是

26、确定的,其余两个未定,所以不能列为第1号。接着检查第2号,由于2号积分器的输入端需要1号提供,故也不能列为第1号。而3号积分器的输入端由2号积分器的初始值便可确定,故列为第1号。这样反复检查,最后得到的计算顺序为345612显然这种计算顺序并非唯一,它和用户编号顺序有关。若把原用户编号作如下改动,即13,24,35,41,52,66则计算顺序为125634虽然顺序不同,只要确定正确无误,那么其计算结果完全相同。3.运行程序块该段主要包括(1运行准备:指针 A 龙格-库塔法系数;B 典型块运行控制。(2典型块计算,目前MCSS仿真程序共有17种典型块,且程序为用户提供了方便的备用模块入口地址,以

27、适应特殊需要。另外程序为了突出简单,积分方法只选择了四阶龙格-库塔法。4.输出数据块MCSS仿真程序追求的是简单,故输出数据的方式很简单,只有三条语句为数据输出语句。每次运行最多可将4个模块的数据输出。在连续时间系统面向结构图仿真中,完全可能出现这种情况,即计算机无法对所有的环节进行排序。换句话说,在系统的仿真结构图中,可能存在一些环节无法安排计算顺序。这种情况我们称系统的仿真结构图出现了代数环。这是需要特别注意的问题,因为这时,系统将不可能进行仿真计算。所谓的代数环就是在系统的仿真结构图中出现了纯粹是有代数环节构成的闭环回路。 环,在该封闭环内没有动态环节,而仅有一些代数环节组成。或者也可以

28、这样说,如果某环节的输出仅仅经过代数运算就反馈到该环节的某个输入,则就构成一个代数环。代数环的解决方法很多,可以根据被仿真对像的具体情况来定。在SimuLink里面,采用了Newton-Raphson技术来求解代数环。由于这种方法利用的是迭代优化,因此存在计算时间长,甚至找不到解的可能性。在此我们介绍两种比较实用的方法。如果系统的仿真结构图比较复杂,并且计算步长也较小,可以在代数环上的任意一个地方增加一个一步延时环。这种方法比较简单,并且有效。 即按上式重新画仿真结构图,就可以避免代数环的出现。因此,一般按标准的一阶微分方程组,即方程右边不出现变量的导数,画出的仿真结构图不会有代数环。本章小结

29、本章是全书的重点之一,所介绍的两种面向结构图仿真方法是目前科学研究和核工程实践中常用的仿真方法。(1当连续系统用结构图形式给定后,离散相似法是一种较为简单的方法,该方法的实质就是在系统必要环节的输入和输出端加入虚拟采样器和保持器,将连续系统离散化,然后分别计算分割开的各个环节的输出量,并按结构图上的关系把相应的输入与输出连接起来,顺序求解计算。由于环节的离散化方程可离线计算,因此该方法突出的特点是运算速度快,但精度低。(2相对于离散相似法而言在本章第四节叙述的连续系统结构图仿真方法应用更广泛。它们的基本原理是一样的,但后者对各环节的计算是在仿真运行中分别计算的,因此它可以提供多种积分解法供用户

30、选择,并且有较高的精度,但带来的问题是计算速度较慢,并且需要排序,以避免结构图中出现的代数环。习题4-1 有一闭环系统如下图所示。 (1求出的,列出求解y(t的差分方程。 (2求出闭环系统的,列出求解y(t的差分方程。(3在闭环入口e处加虚拟采样器及零阶保持器,求出开环的脉冲传递函数W(z,并列出求解y(t的差分方程。(4在系统入口u处加虚拟采样器及零阶保持器,求出开环的脉冲传递函数W(z,并列出求解y(t的差分方程。4-2 已知系统:故差分方 程为试分析步距T应如何选择。若选择的过大,计算时是否会发生不稳定?4-3 试编出下图所示的非线性环节的仿真程序。 4-4 设有系统状态方程如下: (1

31、 (2试编制程序,求出这两组状态方程转移矩阵和。已知采样周期T=0.5。4-5 试用MCSS仿真程序对下图所示系统进行仿真。其中,输入u为单位阶跃函数;参数c1=1,c2=0.05。分析在离散化采样周期(步长T=0.02和T=0.05两种情况下,仿真结果是否相同。 4-6 设有一系统结构如下图所示。当输入为单位阶跃函数时,输出y的精确解为 试用离散仿真程序对上述系统在以下几种条件下进行仿真,并与精确解y(t相比较。(1在离散化采样周期D0=0.02,0.005,0.001,0.002四种情况下的仿真结果与精确解的误差。(2在离散化采样周期D0=0.005条件下,对于把环节合在一起离散化以及两个

32、环节分开离散化两种情况,若将两者所得结果同精确解比较,那一种误差更小一些? 4-7 试用MSCSS程序仿真一个具有延迟环节的系统,其结构图如下图所示。离散化采样周期D0=0.01。 附录:MCSS仿真程序原码500 DATA 1,1,0,0,0510 DATA 2,2,1,8,0520 DATA 3,16,2,0,0530 DATA 4,3,0,3,0540 DATA 5,3,4,5,0550 DATA 6,3,5,6,0560 DATA 7,3,0,6,7570 DATA 8,2,6,7,0580 DATA 0,0,0,0,0590 DATA 1,1,0,0600 DATA 2,1,-1,0

33、610 DATA 3,1,1,0620 DATA 4,0,4,0630 DATA 5,0,-2,0640 DATA 6,0,-10,0650 DATA 7,0,.4,-.1660 DATA 8,1,1,0670 DATA 0,0,0,0680 DATA .05,30,.2,8,2,6,7690 DIM U(200,5,V(60,P(200,3,C(200,E(60,4,Y(200,TT(200,CK4(200,5 700 DIM S(20,200,N(200,ID(20,AD(9,BD(9,CY(9,CU(9710 N9=0:T=0:C(0=0:JNP=1720 PRINT "INP

34、UT I,I6,J,K,L"730 INPUT J,U(J,1,U(J,2,U(J,3,U(J,4740 IF J=0 GOTO 830750 N7=N7+1:I6=U(J,1:IF I6<>4 GOTO 770760 N4=N4+1:ID(N4=J:GOTO 810770 IF I6<>3 THEN 790780 N9=N9+1:U(J,5=N9:V(N9=J:GOTO 810790 IF I6<>1 GOTO 810800 N8=N8+1:N(N8=J:GOTO 820810 U(J,1=-I6820 GOTO 730830 FOR I=1 T

35、O N7 :IF U(I,1>0 THEN 940840 FOR K=2 TO 4:L1=ABS(U(I,K :IF L1<>0 THEN 860850 IF K=4 THEN 930 ELSE 920860 IF N4=0 THEN 880870 FOR J=1 TO N4:IF L1=ID(J THEN 920 ELSE NEXT J880 IF N9<=0 THEN 950890 FOR J=1 TO N9 :IF L1=V(J THEN 920 ELSE NEXT J900 FOR J=1 TO N8 :IF L1=N(J THEN 920 ELSE NEXT

36、J910 GOTO 940920 NEXT K930 N8=N8+1:N(N8=I:U(I,1=-U(I,1:GOTO 830940 NEXT I :N(N7+1=N7+1:GOTO 970950 PRINT "NEXTWORK MUST CONTAIN AT LEAST ONE INTEGRATOR" 960 PRINT "RECHECK DATA STATEMENTS AND RESTART":END970 FOR I=1 TO N7 :IF U(I,1>0 THEN 990980 I4=2:U(I,1=-U(I,1:PRINT "S

37、ORT FAILURE BLOCK"I4990 NEXT I:PRINT "INPUT I,Y(1,Y(2,Y(3"1000 INPUT J,P(J,1,P(J,2,P(J,31010 IF J=0 THEN 10301020 C(J=P(J,1:Y(J=P(J,1:GOTO 10001030 IF N4=0 GOTO 10501040 FOR I=1 TO N4:J=ID(I:P(J,2=I:NEXT I1050 INPUT "D,T2,T1,K1,K2,K3,K4"D,T2,T1,K1,K2,K3,K41060 N1=T1/D:N2=INT

38、(N1+.5:D7=.5*D:D6=D/6:NPRT=INT(T2/T1+.51070 PRINT "TIME"TAB(15;"OUTPUT-"K1;TAB(30;"OUTPUT-"K2; 1080 PRINT TAB(45;"OUTPUT-"K3;TAB(60;"OUTPUT-"K41090 A=11100 B=11110 I7=N(B:I6=U(I7,1:IF I6=0 THEN 16801120 J1=U(I7,2:K=U(I7,3:L=U(I7,4:P1=P(I7,1:P2=P(I7,2

39、:P3=P(I7,31130 C2=C(J1:C3=C(K:C4=C(L1140 ON I6 GOTO 1670,1150,1160,1170,1220,1230,1240,1250,1260,1270, 1280,1300,1310,1420,1470,1520,15701150 C1=P1*C2+P2*C3+P3*C4:GOTO 16601160 I3=U(I7,5:E(I3,A=C2+P2*C3+P3*C4:GOTO 16701170 IF T>D GOTO 11901180 P3=INT(4*P1/D+.5:P(B,3=P31190 IP2=P2:S(IP2,P3=C2:C1=S

40、(IP2,01200 FOR M1=1 TO P3:S(IP2,M1-1=S(IP2,M11210 NEXT M1:GOTO 16601220 C1=SIN(C2 :GOTO 16601230 C1=COS(C2:GOTO 16601240 C1=C2*C3 :GOTO 16601250 C1=-C2 :GOTO 16601260 C1=ATN(C2 :GOTO 16601270 C1=C2/C3 :GOTO 16601280 IF C2<0 THEN C1=C4:GOTO 16601290 C1=C3 :GOTO 16601300 C1=T:GOTO 16601310 IF T>

41、.0001 THEN 13501320 GD5=1:FOR I=1 TO 8:CU(I=0:CY(I=O:NEXT :IF ZD7>.1 THEN 13901330 READ JA,JB:FOR I=1 TO JA:READ AD(I:NEXT I1340 FOR I=1 TO JB:READ BD(I:NEXT I:ZD7=1:GOTO 13901350 IF GD5*P2-T>.000001 GOTO 13701360 IF A=1 GOTO 13801370 C1=CY(1:GOTO 16601380 GD5=GD5+11390 C1=0:FOR I=1 TO 8:C1=C1+AD(I*CU(I-BD(I*CY(I:NEXT I:C1=C1+AD(9*C2 1400 FOR I=0 TO 6:CU(8-I=CU(7-I:CY(8-I=CY(7-I:NEXT I:CU(1=C2:CY(1=C1 1410 GOTO 16601420 IF C2<0 GOTO 14501430 IF C2=0 GOTO 14601440 C1=1 :GOTO 1