控制系统计算机辅助设计与仿真(简

1、华北电力大学自动化系 刘长良 电话Email: changliang_liu,控制系统计算机辅助设计与仿真,教材,1、韩璞，自动控制系统数字仿真，中国电力出版社，1996 韩璞，控制系统数字仿真技术，中国电力出版社，2007 2、薛定宇，控制系统仿真与计算机辅助设计，机械工业出版社，2005,参考教材,1、薛定宇，控制系统计算机辅助设计，清华大学出版社，1996 2、孙增圻，控制系统计算机辅助设计与仿真，清华大学出版社，1988 3、刘金琨，先进PID控制及其MATLAB仿真，电子工业出版社，2003 4、张葛祥，MATLAB 仿真技术与应用，清华大学出版社，2003

2、,一、本课程的特点及要求,本课程是一门综合应用课程。,目的：在自动控制原理、过程控制的基础上，进一步掌握控制系统计算机辅助设计与仿真的有关理论与方法以及应用这些理论与方法进行控制系统研究的实用技能。 掌握：控制系统数字仿真、控制系统计算机辅助设计、以MATLAB语言为基础的控制系统计算机辅助设计与仿真的程序设计方法以及这些方法在工程实际中的应用。 培养：运用相关理论和计算机辅助手段解决实际问题的能力。,第一章 绪论,计算机辅助建模 计算机仿真 控制系统计算机辅助分析 控制系统计算机辅助设计,二、主要内容,三、动态过程数学模型概述,1. 数学模型: 是对于现实世界的一个特定对象，根据其内在规律，

3、做出一些必要的假设，运用适当的数学工具，得到一个数学结构。 数学模型又可分为静态模型与动态模型 静态模型：用来描述系统在稳定状态或平衡状态下各种输入变量与输出变量之间的关系。 例如：当机组运行在稳定状态时，输入的物质及能量保持不变，机组各系统的参数也将保持稳定，这些稳定工况下各参数之间的关系便可用静态模型描述。静态模型主要用于机组的设计计算及校核计算，一般要求具有较高的精度。,动态模型： 用来描述系统在过渡过程中各种变量随时间变化的关系。 当系统从一个稳定状态变化到另一稳定状态时，哪些参数会发生变化，其变化的速度及历程如何，这些都属于动态模型要研究的问题。 例如，当燃料量变化时，机组原来的平衡

4、状态就会受到破坏，电功率等参数都将发生变化，经过一段时间运行，机组又将达到新的平衡状态。这个动态过程中电功率的变化规律需要用动态模型描述。,2、 动态过程建模的一般方法 (1). 机理法白盒法 依据基本的物理定律 优点：范围广、物理意义明确 缺点：精度差，与简化程度有关，与过程的不确定性有关。 (2). 试验建模黑盒法 系统辨识 时域法： 飞升曲线法 频域法： 正弦信号、离散信号 相关分析法：随机信号 参数估计法：最小二乘法、极大似然法等 优点：与实际设备特性一致，模型精度高，包含了各种干扰模型 缺点：需要大量的人力物力，模型的应用应限制在试验范围内。 (3). 灰盒法 ： 先建立机理模型，再

5、进行辨识、修正模型结构和数据,3. 火电厂热工过程建模与仿真,火电厂热工过程特点： 系统复杂、多变量、大迟延、非线性、时变性 建立模型时要有针对性： (1) 要明确建模的目的及模型的用途； (2) 确定建模的方法 (3) 对系统进行合理的简化 (4) 确定模型的输入、输出 (5) 建立模型 (6) 分析、验证、修正模型 (7) 模型应用,火电厂设备,计算机接口,实际DCS操作,数学模型,计算机接口、网络,仿真DCS操作,实际电厂,仿真机,机理法建模实例：弹性阻尼系统建模实例 问题：研究行使在不平路面上的汽车的颠簸情况。 方法：物理模型研究 过程：（1）采用类比方法得到相应的概念模型 （2）建立

6、数学模型 建立弹性阻尼系统模型过程： （i）分析 （ii）假定 K、B线性 （iii）建模：牛顿第二定律 得到弹性阻尼系统模型：,四、数字仿真的基本概念 数字仿真：求解数学模型数值解的过程 仿真分类：按模型类型分；按系统的连续性分；按时间尺度分； 实时仿真、非实时仿真,弹性阻尼系统建模实例,（3）模型求解,（i）模型变换 转换为状态方程：,初始条件：x1(0)=0, x2(0)=0;,（ii）选择数值计算方法得到仿真模型 采用欧拉法可写出迭代公式如下：,iii）选择编程语言，编程并进行仿真研究 例如：Matlab、C、Basic语言等 （4）对仿真结果进行分析并验证模型 一般需要试验数据的验证

7、，或采用实验数据修正模型 （5）模型应用，在确认模型正确后，可应用该模型对弹性阻尼系统进行研究,物理原型概念模型（分析、类比得到） 数学模型（假设，采用物理定律得到） 仿真模型（采用一定的数值计算方法） 数字模型（程序代码、通过编程得到） 实用模型（多次的模型分析、检验、修 正） 模型应用（通过模型研究物理系统的特性） 建模与仿真需要明确的问题： 明确建模目的 模型的输入、输出 模型的描述形式 建模对象机理及特性分析 编程语言的选择 模型的验证方法：静态工况、动态过程,动态过程建模与仿真研究的一般步骤,仿真的应用领域，发展史，火电站仿真 机的发展及现状 航空航天 核电、火电 交通、运输 通讯、

8、电子 经济、人类 我校火电仿真技术的发展 1111111111111,五、 计算机辅助分析 、数学模型之间的相互转化 、控制系统稳定性分析 频域法：由W(S)分析系统的稳定性， 劳斯判据、波特图、奈奎斯特图、 根轨迹等 时域法：状态空间法。,六、 计算机辅助设计 频域法： 基于传递函数的控制器设计方法； 例如：串联校正器的概念及设计方法； 极点配置设计方法； 时域法： 基于状态空间模型的控制器设计方法； 线性二次型最优调节器的设计方法； 观测器的基本设计方法； 过程控制系统：PID控制器设计及参数整定方法；,机理法建模： 实例4、椅子着地问题： 4条腿长度相等的椅子，放在起伏不平的地面上，4条

9、腿是否能同时着地？,假定地面为连续平面，则转动椅 子在90内，总能找到一点使四 条腿同时着地。 让椅子原地旋转，x为旋转角度。 设f(x)为A、C两腿距地面距离 之和，g(x)为B、D两腿距地面距 离之和。,由于任意三腿总在一个平面上，所以任意地点三腿可同时着地。总有f(x) , g(x)之一为0，即f(x)g(x)=0 数学描述： 已知：f(x) , g(x)为x的连续函数，已知g(0)=0,f(0)0 且 f(x)g(x)=0 求证：存在x0 使得g(x0)=0 证明：令h(x)=f(x) - g(x), 则h(0)=f(0) - g(0)=f(0)0 将椅子转动90度，对角互换，则：由g

10、(x)=0,f(x)0 可得：g(90)=f(0)0, f(90)=g(0)=0 即：h(90)= f(90)-g(90)状态方程 方法（1）：传递函数 微分方程状态方程,例3：求下列传递函数对象的状态方程描述,解：先将传递函数进行拉氏反变换得到微分方程,再将微分方程转换成状态方程,例4：,2. 传递函数状态方程 方法（2）：直接转化 串联法,等效框图：,解：,同理：,状态方程：,输出方程：,初值：x1(0)=0; x2(0)=0; x3(0)=0;,2. 传递函数状态方程 方法（3）：直接转化 并联法,等效框图：,并联法的特点：方法通用，A为对角阵,同理：,状态方程：,输出方程：,初值：x1

11、(0)=0; x2(0)=0; x3(0)=0;,注意： 1. 状态方程描述不是唯一的，不同变化方法得到的状态方程及输出方程是不一样的 ，但阶次不变。 2. 要想得到A为对角阵，必须采用并联法！,例如：前述的例题也可以采用并联法求解,3. 传递函数框图形式的描述状态方程,例5：求如下控制系统的状态方程描述,等效框图：,按照前述方法，可得： x1=0.04e， x2=-0.1x2+0.1e1 x3=-0.05x3+0.05x2 e=u-x3， e1=0.05e+x1=0.05u-0.05x3+x1,把e,e1代入上面各式，整理可得状态方程描述:,状态方程：,输出方程： y=x3,初值：x1(0)

12、=0; x2(0)=0; x3(0)=0;,本章复习题,2. 将微分方程转为状态方程描述，假定初值为0,3. 将传递函数转为状态方程描述，要求A为对角阵。,简述热工过程机理建模的基本假设，并给出常用的公式,4. 求以下控制系统状态方程描述，输入为u,输出为y,假定初值为0.,第三章：数字仿真算法-数值积分法,问题： 已知状态方程，如何求状态变量的数值解？,已知： 状态方程： X=AX+Bu 输出方程： y=Cx+Du 初值： X(0) 求：各个时刻输出y的值 即：给定一个时间计算步距DT，求y(DT), y(2DT), y(3DT), y(4DT),.y(kDT).y(nDT),例：弹性阻尼系

13、统数值计算,设状态变量：x1=x, x2=x 得：,X1=x2, x2=-k/m x1 b/m X2 +1/m u,如何计算x1,x2在各个时刻的数值？,对于：,在某个时刻t0, 近似有：,由此可得：,即：对于微小的DT，可以由t0时刻的x1,x2, 近似计算t0+DT时刻的x1,由前面，已知x1(0) ,x2(0), 这样有如下递推公式：,.,对于任意一点， 有如下递推公式：,简记为：,按此递推公式即可求得各个时刻x1的数值,一、数值积分法仿真算法,1. 数值积分仿真算法原理,对于： X=AX+Bu 设 F(t)=A X + Bu，则：X=F(t),两边积分,得X计算通式：,取t=kDT,

14、t=(k+1)DT,由高数可知： 在kT,(k+1)DT区间，存在一点，F*,可得：,如何计算积分项是数值积分法的关键问题！,如何计算F*是关键问题！,2. 各种数值积分仿真算法,（1）欧拉公式：,取 F*=F(kDT) , 得：,简记：,注意到：,F(t)=A X + Bu,得：,由此，得到欧拉公式仿真算法：,问题： 在区间KDT，(k+1)DT内F(t)是变化的，近似取F(t)=F(KDT)会造成较大的误差。,欧拉公式几何解释：,精确解：,欧拉公式： 用矩形面积近似代替积分面积,DT*F(K),误差,分析：欧拉公式简单实用，但当DT较大时，会造成较大的误差。要想提高精度，必须选择较小的计算

15、那步距DT。DT越小，精度越高。,DT的提高受计算机精度、速度的限制，在DT一定的情况下，如何提高计算精度？,问题：,（2）梯形公式：,用梯形面积近似代替积分面积,取：,由此，得到梯形公式仿真算法：,梯形公式存在的问题：,计算F(k+1)时需要用到X(k+1),无法使用！,怎么办？,放弃,改进,如何改进？,把欧拉公式引入梯形公式，先采用欧拉公式预测F(k+1),再采用梯形公式计算。,改进方法：,改进后的梯形公式：,预报校正公式,改进的梯形公式几何解释,在DT一定的情况下，能否进一步提高计算精度？,问题：,欧拉公式：只采用k点的信息，计算F*,算法分析：,梯形公式：采用了k点和k+1点的信息，计

16、算F*,是否可以采用k点，k+1,k+1/2点的信息计算F*,进一步提高精度？,B、采用欧拉公式预测k+1点的F值F(k+1),C、计算F*:,D、采用欧拉公式计算，X(K+1)=X(k)+DT F*,梯形公式计算步骤：,A、先计算k点的F值F(k),（3）四阶龙格库塔法：,2. 根据F1预测k+1/2点的F值F(k+1/2), 记做F2,6. 根据F*, 采用欧拉公式计算，X(K+1)=X(k)+DT F*,计算步骤：,1、先计算k点的F值F(k)，记做F1,3. 根据F2重新预测k+1/2点的F值F(k+1/2), 记做F3,5、计算F*:,4. 根据F3预测k+1点的F值F(k+1),

17、记做F4,计算步骤：,1.,2.,3.,4.,5.,6.,总体截断误差： 1、欧拉公式 O(DT) 2、梯形公式 O(DT2) 3、四阶龙格-库塔法 O(DT4),3. 仿真算法精度分析：,算法精度：,总体截断误差,局部截断误差,4. 仿真算法稳定性分析,考虑较为简单的情况，对于一阶系统,若采用欧拉公式计算，公式如下：,由k时刻误差引起的k+1时刻的误差计算如下：,若|1+DTa|1，注意一般a0，例如：1+DTa = -2，则：,由-1-2, a0, 得收敛区间：DT-2/a,反之，若|1+DTa|1 则数值计算收敛,显然，数值计算发散,精度与稳定性,需要注意： 1.精度与稳定性是两个不同的

18、概念,两者之间没有 直接的联系； 2.上述每种数值积分算法都有自己的收敛区间； 3.收敛区间与算法有关，也与对象特性有关; 4.仿真步距（计算步距）DT，越小，算法稳定性越好；,5. 其它仿真算法,（1）三阶龙格-库塔法,以上数值计算方法为单步法，由 K 时刻的值推导K+1 时刻的值。,（2） 阿达姆斯法,多步法：,由 K，K-1，k-2.时刻的值推导K+1 时刻的值。,误差 :O(T4),单步法:,阿达姆斯法：,（3）隐式欧拉公式,若取 F*=F(k+1) , 得：,对于高阶系统需要求逆矩阵，采用高级语言较为复杂，但对于单个微分方程：,把F(k+1) 代入得：,隐式欧拉公式如何求解？,隐式欧

19、拉公式的优点：,稳定性好,对仿真步距DT没有要求，不管DT取多大，计算不会发散,注意：虽然不会发散，但误差可能很大！,说明：精度与稳定性是两个不同的概念,二、仿真程序设计,例：对于前述弹性阻尼系统模型,采用高级语言编程仿真，例如Matlab，C，Basic等,采用欧拉公式如下：,#include ； #include ； #include ； float K=0.7,B=0.8; X1 = 0, X2 = 0, M = 1, u = 1; int i,NP; float y1000,F1,F2, ST=20,DT=0.1; main() NP = (int)(ST / DT);,C语言仿真程序

20、设计,for(i=0;ihelp AB 显示： 弹性阻尼系统仿真程序 函数格式调用为ABC(ST,DT),（8）绘图函数,x=1,2,3,4,5,6; y=110,120,150,120,145,133; plot(x,y),注意： x, y数据个数必须相同,Plot( (1:NP)*DT,y) 等价于,T= (1:NP)*DT; Plot( T,y),for i=1:NP F1=X2; F2=-K/M*X1-B/M*X2+u/M; T(i)=i*DT, Y(i)=X1; X1 = X1 + DT * F1; X2 = X2 + DT * F2 end Q=x1; plot(T,Y); tit

21、le(弹性阻尼系统);,程序,（9）其它绘图函数,grid 网格,title(图形的名称);,hold on 锁住屏幕,xlabel (时间 t/s),ylabel (位移 y/cm),四、热工过程仿真,质量守恒应用1：求水箱水位,A. 假定：,B. 模型输入：u1，输出：L,C. 动态模型：,D. 仿真模型：,显式欧拉公式,研究重点：系统从一个稳定状态变化到另一个稳定状态时所经历的历程。,E. 仿真研究,（1）确定系统初始稳定状态； （2）改变系统的输入，观察输出的响应过程,假定初始稳态如下：L0=100 cm, F=900 cm2 , W10=W20=900 g/s, ro=1 g/cm3

22、, u10=0.5 。,由此可计算得到动态模型中的系数： K1=1800, K2=90。,仿真时假定阀门开度由50%变到60%，相应的matlab仿真程序如下：,clear; L=100; F=900;K1=1800; K2=90; RO=1; u1=0.6; ST = 2000; DT = 10;NP = ST / DT; for i=1:NP W1=K1*u1; W2=K2*sqrt(L); Y(i)=L; L=L+DT*(W1-W2)/(RO*F); end plot(1:NP)*DT,Y); title(单容水箱仿真);,水箱水位仿真研究,水箱水位仿真曲线,DT仿真研究-显式欧拉公式,

23、水箱水位仿真研究-隐式欧拉公式,D. 仿真模型：隐式欧拉公式,非线性方程，无法直接求得L(k+1)的表达式,需要把上式中的非线性部分进行线性化后，才能得到L(k+1)的表达式。,常用处理方法：,即：把原来非线性项分为两部分：系数部分和线性部分。系数部分用K时刻的值近似计算，得到B2(k),整理得：,for i=1:NP B2=K2/sqrt(L); L=( RO*F/DT*L+K1*u1)/(RO*F/DT+B2); Y(i)=L; end,不同仿真步距是输出曲线： DT=10s, 500s,1000s,2000s,5000s,隐式欧拉公式仿真曲线,该系统由两个水箱组成，是二阶系统，所以其动态

24、过程应包括两个动态方程。,对于水箱1、2分别列出其动态模型如下：,多容水箱系统,例：如图水箱系统：,显式欧拉公式：,仿真模型：,仿真过程如下：,（1）确定初始平衡状态； 假定L10=100 cm, F1=900 cm2 , W10=W20=W30=900 g/s, Ro=1 g/cm3, u10=0.5, L20=64 cm, F2=600 cm2 （2）计算动态模型中的系统常数： 根据上面的初始条件，计算可得：K1=1800, K2=150, K3=112.5。 （3）假定调节阀开度u1由0.5变成0.6，研究动态过程,不同仿真步距是输出曲线：显式欧拉,隐式欧拉公式仿真模型如下：,按照同样的

25、线性化方法可得：,整理得仿真模型：,不同仿真步距是输出曲线：隐式欧拉,如图所示：流体网络中一个压力节点,建模时假定： (1) P点周围的管道容积都集中在P处 (2) 忽略流体密度变化,根据质量守衡可得：,二、 质量守恒应用之二：流体网络压力节点模型,对于热工过程而言，一般情况下，温度的变化速度要比压力的变化速度慢得多，忽略温度的变化,K1的大小反映了流体的可压缩性，K1值越大，说明流体可压性越强，反之，可压性则越小，对于不可压流体K10。 需要注意的是，只有理想流体才可能是不可压流体，实际流体都是可压的。,忽略流体密度变化，则各路流量与压差关系如下：,若采用显式欧拉公式：,若采用隐式欧拉公式：

26、,按前述方法对上式进行线性化，设：,压力节点仿真模型：,考虑流体网络中各节点的计算顺序，在计算P(k+1)时，P1(k+1)、P2(k+1)、P3(k+1)可能还是未知，为使模型通用，所有相临节点压力都取为K时刻的值，即：,如图， 对于小流量，为简单起见，直接引用流量,若采用隐式欧拉公式：,对于任意多支路的压力节点，假定周围有i个压力节点，j路微小流量：,建模时假设： （1）外壁绝热； （2）采用集总参数法，假定金属管道壁温均匀一致； （3）忽略工质定压比热Cp随温度的变化,三、 能量守恒应用之一：绝热蒸汽管道模型,根据质量守衡：,根据能量守衡,对于工质有：,对于金属有：,工质和金属之间的换热

27、量：,对流换热系数：,1. 基于温度的集总参数模型,由于实际过程中压力、流量变化较快，而温度、焓值变化较慢，建立模型时一般把流动过程与传热过程分开考虑。质量守衡动态过程在压力节点模型或水箱水位模型中考虑，能量守衡方程在温度、焓值模型中考虑。,取蒸汽侧集总参数T=Tout ，采用显式欧拉公式得：,若采用隐式的欧拉公式，注意到Q与T、Tm有关，为了增加模型的稳定性，需要将Q带入,考虑到计算顺序，前面公式中用到的Tm(k+1)要用Tm(k)代替才能计算,由于水蒸汽的定压比热随压力、温度的变化很大，定值假定会造成较大误差，在建模时常采用焓值作为热平衡计算的基础，建立基于焓值的集总参数模型:,H为集总参

28、数，取H=hout。由于计算Q时要用到蒸汽温度，还需要根据焓值计算温度，然后再计算传热量Q。,假定蒸汽压力为P,则：,2. 基于焓值的集总参数模型,由于Q与H没有直接的函数关系，所以无论用显式还是隐式欧拉公式，只能采用Q(k),Tm的递推公式同前：,在该模型中，若考虑由于外壁与环境的自然对流换热，则可假设总的散热系数为Kamb，环境温度为Tamb，模型中蒸汽侧不变，只需在金属侧减去管道对环境散热量 ：,计算Tm的递推公式也应进行相应的改动。,建模时假定：,（1）整个换热器用一段逆流管代替，沿管长方向吸热均匀，金属温度均匀一致，为Tm ； （2）忽略烟气定压比热随温度的变化； （3）忽略流动的动

29、态过程，取Ws2=Ws1=Ws, Wg2=Wg1=Wg ，工质压力为Ps；,四. 能量守恒应用之二：单项介质换热器模型,烟气侧：,动态模型如下：,工质侧：,金属侧：,换热计算中蒸汽、烟气取平均温度，采用隐式欧拉公式：,整理得：,输入：Ws, hs1, Ps, Wg, Tg1 系数：Kg1, Kg2, Ks1, Ks2, Mg ,Cpg, Mm ,Cm, Ms 输出：hs2, Ts2, Tg2, Tm,整理后模型的一般形式,整理好的仿真模型,整理好的仿真模型如下,混合联箱：联箱内有一电加热器，热功率为Q,建模时假定： （1）忽略流体密度变化，水箱等截面，面积为F； （2）集总参数法，整个联箱内工

30、质温度均匀一致； （3）忽略工质定压比热Cp随温度的变化；,五、 质量守恒-能量守恒同时应用：混合联箱模型,质量守恒：,能量守恒：,虽然上式中V与L有关，但建立动态模型时一般并不将上式按微分式子展开，而是处理（2）式先将L按常数处理，得到仿真模型，再和（1）一起迭代运算。,动态模型：,微分项的近似处理方法：,简化如下：,集总参数T取为出口温度t3，采用隐式欧拉公式：,仍按前面的线性化处理方法可得隐式欧拉公式模型:,模型输入: W1，t1, W2, t2, Q,模型输出: L，t3,已知：,第四章、离散相似法,一、离散相似法原理,状态方程：,输出方程：,初值：,求：Y的数值解,方法一：前述的数值

31、积分法,方法二：离散相似法,通解：,K时刻：,K+1：,如何计算eAT及积分是该算法的关键！,递推关系：,与数值积分法比较,数值积分法：近似处理F,离散相似法：近似处理,取t1=(K+1)T-t,进行积分变换：,零阶保持器离散相似算法如下：,在区间kT,(k+1)T, 近似取U(t)=U(kT),零阶保持器,（1）、计算(T) ：,一般指数 ：,矩阵指数 ：,近似取前几项,（2）、计算m(T) ：,A、当系统为一阶系统时，A变为常数，求(T)、 m(T)变得较简单； B、 MATLAB有相应的函数,(T)、 m(T)的计算很复杂！,在区间kT,(k+1)T,近似取,代入得：,整理得：,其中：,

32、三角阶保持器,三角保持器离散相似算法如下：,方法三、通过拉氏变化求解：,求拉氏变换： sX(s)-X(0)=AX(s)+BU(s) (sI-A)X(s)=X(0)+BU(s) X(s)=(sI-A)-1X(0)+ (sI-A)-1BU(s) 求拉氏反变换 X(t)=L-1(sI-A)-1 X(0)+L-1(sI-A)-1BU(s),根据卷积定理：,对于：,设：,由此可得：,与离散相似法比较：,可得：,即：可以通过拉氏反变换计算,同理：x2=-1/5x2+1/5x1 e=u-x2 ,整理得：,解: (1) 先求系统的状态方程 设定状态变量x1,x2,例1：对如下系统进行仿真,x1=-1/4x1+

33、1/4e,求(T),m(T), n(T),由于不是对角阵，求逆矩阵比较困难!,若A是对角阵，例如,可见，离散相似法过于繁琐!,解决方法1：,系统为一阶系统时，A=a, eAT=eaT,对每个一阶系统分别进行离散相似,x1=-1/4x1+1/4e,A=-1/4, B=1/4, (t)=e-t/4, (T)=e-T/4,x1(K+1)= e-T/4x1(K)+ (1-e-T/4)e(K) x2(K+1)= e-T/5x2(K)+ (1-e-T/5)x1(K),同理：,e(K)=u(K)-x2(K) x1(K+1)= e-T/4x1(K)+ (1-e-T/4)e(K) x2(K+1)= e-T/5x

34、2(K)+ (1-e-T/5)x1(K) y(K+1)= x2(K+1),仿真模型如下：,（1）原来是对u进行近似处理，这次分别对两个子系统的输入e,x1进行了近似处理，误差增大。 （2）当u是阶跃输入时，对u进行近似处理不会引入误差，而第二种方法对e,x1近似处理，引入了两次误差。,注意：,分析发现，计算x1(k+1)用到e(K+1) e(K+1)=u(K+1)-x2(K+1) 可见，该模型是无法使用的。,为了减小误差，采用三角保持器：,即：,A在区间kT,(k+1)T,整理得：,处理方法,方法1：对系统1（与反馈环节相邻）使用零阶保持器,方法2：改进三角保持器,一阶保持器,B在区间kT,(

35、k+1)T,取,这样近似处理后，算法实用化。,滞后一拍三角保持器,整理得：,二、保持器的概念,b. 一阶保持器,c. 三角保持器,d. 滞后一拍的三角保持器,采样系统：采样开关、保持器 采样处理离散，信号恢复（保持）相似,2保持器种类及频率特性,1. 离散相似概念的来源：,a. 零阶保持器,3. 各种离散相似处理后信号的变化,（1）.多次离散后造成信号有较大的相位滞后,（2）.多次离散后造成信号有一定的幅值衰减,为了提高精度，对离散相似信号补偿，包括对幅值及相位的补偿,为了简单,一般取,4. 补偿器,补偿器：,即:输入信号超前一拍,例如对于零阶保持器,加入补偿环节后:,小结：,1. 对u离散相

36、似，各种算法都能使用，但需要求得整个系统的状态方程； 2对于负反馈系统，对e离散相似，三角保持器不能使用，对x1离散相似，各种算法都能使用； 3. 对e离散相似, 不能使用补偿器； 4. 零阶保持器可以无失真的再现阶跃信号； 5. 其它保持器可以无失真的再现斜波信号； 6. 对各个子系统分别进行离散相似可以简化仿真算法，但会引入新的误差；,离散相似法程序设计,采用matlab 函数进行弹性阻尼系统仿真的 程序如下(离散相似法零阶保持器) clear all；clc； m=1;b=0.5;k=1;u=1;X=0 0;T=0.01;ST=30; A=0 1;-k/m -b/m;B=0 1/m;NP

37、=round(ST/T); Q=expm(A*T); E=inv(symsub(symmul(s,eye(2,2),A); F=ilaplace(E) G=symmul(F,B),Qm=numeric(int(G,t,0,0.01) for i=1:NP y(i)=X(1); X=Q*X+Qm*u; end plot(0:NP-1)*T,y,y);,对于单容水箱系统：,需要注意，上述算法是对于线性系统推导得出的，只能用于线性系统。,A=-B2/(roF)， B=K1/(roF) (T)=eAT,m(T)=K1/B2(1-eAT) L(K+1)= (T)L(K)+ m(T)u(K),clear;

38、 L = 100; F=900;K1=1800; K2=90; RO=1;u1=0.6; ST = 2000; DT =100; NP = round(ST / DT); for i=1:NP W1=K1*u1; B2=K2/sqrt(L); AA=-B2/(RO*F); BB=K1/(RO*F); FA=exp(AA*DT); FAM=K1/B2*(1-FA); Y(i)=L; L=FA*L+FAM*u1; end plot(1:NP)*DT,Y,r); title(采用离散相似法进行单容水箱仿真);,仿真结果分析： 离散相似法与 隐式欧拉公式对比,三、Z变换法,原理：脉冲传递函数直接对应差

39、分方程,1. 脉冲传递函数的求取：输入输出的z变换之比,拉普拉斯变换定义:,：拉普拉斯算子，,z变换定义：,已知X(t)的采样信号x*(t),对于kQBQC, 即QA最差，QB为除QA之外的最差点。 则新点A1: x10*=x11+x12-x10; x20*=x21+x22-x20; if (QA-QC)/QAEPSQ, 结束。,(3) 单纯形的压缩： A、若QA1QB, 则新三角形中A1最差，若按前面方法求反射点，则反射点为原来的z1。此时，说明新点沿反射方向前进太多，应适当减少前进距离。 即取A2: x10*-x1c=k(x1c-x10) (-1k1 ) x10*=(x11+x12)*(k+1)/2-kx10; 同理：x20*=(x21+x22)*(k+1)/2-kx20;,4) 在(3)B 中取k=0.5若QA2QB,取k= - 0.5 求得A3点的QA3重新判断。,(5) 若QA3QB,取k=0 求得A4点, 无法形成三角形 可构造新点C1，x1c1=x12+delta1*lx1, x2c1=x22+delta1*lx2, 由C1,B,C三点构成新三角形，转入（2）。,编程改进： 1。为保证三角形的性状，将dta,Ti寻优区间等量化，例如：dta=0.1,10,ti=1,10