火电厂仿真(课件)

火电厂仿真太原理工大学电气与动力工程学院热能系1参考书：

张家琛编《火电厂仿真》

北京：水利电力出版社1994、10

唐世林编《电站计算机仿真技术》

北京：科学出版社1996、11

2前言

仿真：模拟真实设备运行模拟真实的过程模拟真实的环境

3一、火电厂仿真技术应用的优越性

1.减少培训费用，缩短培训周期仿真装置可模拟现场机组的启停、运行操作，还可进行现场不能进行的事故演练和操作实验、分析事故及电厂运行的经济性，不仅可对新的运行人员进行教育，而且可对原有人员进行再教育，对提高工作人员的素质，提高发电站运行水平和管理水平有重要意义。42.节约研究开发时间和经费利用仿真装置可不受现场各种条件的约束，对改造旧的和开发新的动力装置有重要意义。3.提高机组的安全性、效率和使用寿命由于培训和研究开发试验不在现场进行，可以避免现场的人为失误、减少停机次数、提高设备的安全性、提高机组使用效率及寿命。51955年，Profon首先将仿真技术应用于锅炉研究，其研究分为建模理论、燃烧模型、锅炉部件模型三部分，相继发表了许多有意义的研究成果。——可谓电站仿真的开端。1958年，K·L·Chien建立了第一个汽包锅炉模型。随着电力工业的发展和技术的进步，发电机组的容量不断扩大、参数不断提高、自动化程度越来越高，对运行操作人员的技术知识、操作水平、应变能力、熟练程度提出了更高的要求。为此，从20世纪60年代中期，研究人员开始探索采用脱离实际发电机组的实时仿真装置，对运行人员进行培训。二、火电厂仿真的发展历史61968年，美国建成的Dresden电站仿真装置，标志着电厂仿真应用的开始。在1970年以前，限于数字计算机的水平，当时开发的仿真装置都是采用的模拟电路，且只能对电厂局部范围进行仿真。进入20世纪70年代，水蒸汽状态方程标准的制订、发电机组数学模型的完善与优化，数字计算机技术的进步，促进了火电厂仿真的深入。1971年，美国Singer公司Link仿真系统分部、EAI公司和英、日等国几乎同时推出采用数字计算机或数字—模拟混合机的火电机组全仿真系统。在培训中收到良好效果。7随后，各国一些大的电气或电子公司相继成立了仿真公司或分部，专门研究和开发仿真系统。例如：美国：Singer公司（SingerCo·LinkSimulationSystemsDivision）EAI公司（ElectronicAssociates，Inc）CE公司（CombustionEngineering）EPRI公司（ElectricPowerResearchInstitute）Gould公司、Ominidata公司、Autodynamics公司加拿大:CAE公司（CanadianAviationElectronics）英国：CEGB公司（CentralElectricityGeneratingBoard）Link—Miles公司德国：KruppAtlasElectronic芬兰：NokiaElectrnics我国:

清华大学1974年开始了火电机组仿真系统的研制工作，

81975年，美国制订了全复制型核电站培训仿真器国家标准（Nuclearpowerplantsimulatorsforuseinoperationtraining）1979年3月28日美国三里岛核电站（ThreemilesIslandＮpp）由于运转员的失误发生了严重事故。调查结果表明，在核电站操作员的培训方面存在着严重问题。1980年3月美国NRC(NuclearRegulatoryCommission)发布了一套考核和鉴定反应堆操作员的新标准，它规定了反应堆操作员的起码资格。91981年8月7日美国NRC规定：从即日起将把仿真器作为发放运行许可证审批的一个组成部分。同时美国能源局规定：①未取得仿真培训合格证的任何人员不许上岗进行核电站任何操作。②核电站实时仿真培训系统应尽量是“完全复制控制室”的全仿真系统。101983年5月，清华大学开发的我国第一台200MW燃煤机组部分范围仿真机通过鉴定，并投入使用。1985年，美国经两次修订，公布了新的全复制型核电站培训仿真器国家标准（ANSI/ANS3·5，1985）1986年，华北电力学院从美国引进的550MW部分范围仿真机投入使用。1988年1月，从美国SingerLink公司引进的第一台300MW燃煤电厂全范围培训仿真机投入使用。111988年，清华大学开发的我国第一台200MW火电机组全范围培训仿真机通过验收，并投入使用。南京工学院开发完成125MW火电厂仿真机。1989年，华电计算机仿真公司、华北电力学院仿真与控制技术研究所、西安热工所先后开始开发300MW火电厂全范围培训仿真机。西安热工所还着手开发500MW火电厂仿真机。1982年12月，美国颁布新的《核电站操作员培训仿真机》标准ANSI/ANS3·51990年8月，美国颁布《火电厂仿真机功能要求》ISA-77·20标准。121991年，华北电力学院仿真与控制技术研究所开发的我国第一台300MW燃煤机组全范围培训仿真机在银川投入使用。1993年，我国能源部提出了《大型火电机组仿真培训装置技术规范》（试用稿）。1993年4月，华电计算机仿真公司开发的300MW机组全范围仿真机在武汉电力学校投入运行。1993年6月，美国又通过了《火电厂仿真机功能要求》ISA-77·20-1993版本标准。1994年，西安热工所开发的我国第一台500MW火电机组仿真机在太原电力高等专科学校通过验收。131995年8月，由华北电力学院开发的我国第一台500MW超临界火电机组全范围培训仿真机在华北电管局盘山电厂投入使用。1995年底，我国已投入运行和正在制造的200MW以上火电和核电机组培训仿真机达到60多台；14第一章计算机仿真概述第一节研究问题的方法一、种类：直接试验物理模型试验（一次模化技术）数学模型计算机仿真模型试验（二次模化技术）15二、比较直接试验：耗资大、周期长、安全性差。一次模化法：耗资较大，不在现场、安全性好，具有一定的灵活性。计算机仿真：适应性强，灵活性大（改变软件或参数可对不同问题进行拟）安全性好，准确性待验证。161、物理仿真：按真实系统的物理性质构造系统的物理模型，再现系统的一些特性。例如:新型锅炉生产前，各部件按一定比例缩小，构成一个几何相似的物理模型。对此分析研究可发现该锅炉在结构和形态上是否合理，安装、操作和检修是否便利。这种模型仅反映了其结构特性，而未能反映锅炉内部传热学、热力学、流体力学的特性。一、计算机仿真的概念第二节计算机仿真技术简介172、数学仿真：按真实系统的数学关系构造系统的数学模型，即将实际系统的运动规律用数学形式表达出来，并在数学模型上进行试验，再现系统的某些特性。数学模型能精确反映系统内部的各种静态和动态特性，如锅炉运行中的燃料化学反应、传热过程、能量储存与释放、工质循环流动的特性。数学仿真一般是利用计算机对系统数学模型进行运算和试验的。因此，利用计算机实现的系统数学仿真也称为：计算机仿真。183、物理-数学混合仿真：将真实系统的一部分用数学模型描述，并放到计算机上运行，而另一部分则构造其物理模型或直接采用实物，然后将它们连接成系统，并在此系统上进行试验，再现系统的某些特性。这种仿真又称：半物理仿真。194.计算机仿真定义：将一个能近似描述实际系统的数学模型进行第二次模化，转变为一个仿真模型，并将其放到计算机上进行运行，以再现系统某些特性的过程。第一次模化：研究实际系统与数学模型之间的关系；第二次模化：研究数学模型与计算机之间的关系。计算机仿真是本课程讨论的主题。205.计算机仿真系统——用于系统仿真的一套软、硬设备。——其主体是计算机。按计算机的类型分，有以下仿真系统：1、模拟仿真系统2、数字仿真系统3、数模混合仿真系统21模拟仿真系统——用模拟计算机组成的仿真系统。优点：运算是“并行”的、“连续”的。因此，运算速度快，且更接近连续系统。缺点：运算精度低、线路复杂、对采样和逻辑系统的仿真比较困难、仿真的自动化程度低（依靠排题板接线）。22数字仿真系统——用数字计算机组成的仿真系统。优点：被仿真的系统包含在一组程序中，仿真自动化程度高，使用方便；运算精度高；易实现逻辑处理和非线性环节；程序和参数修改容易。缺点：运算过程是“串行”的，运算速度相对较低、实时仿真和寻优计算等不如模拟仿真系统快。用微型数字计算机阵列组成的仿真系统，称为全数字仿真系统23数模混合仿真系统——用模拟计算机和数字计算机组成的仿真系统。由于模拟计算机和数字计算机的优缺点是互补的。因此，该系统达到了扬长避短的目的。该系统适用于：（1）要求与实物连接，又有许多复杂函数需计算的实时仿真；（2）需要进行反复迭代计算（如统计分析、参数寻优）的仿真；（3）对计算机控制系统的仿真；24二、计算机仿真的基本步骤（1）根据系统研究的目的和范围，建立系统的数学模型；（2）将数学模型转换成计算机能接受、可运行的仿真模型；（3）根据仿真模型编写仿真程序；（4）对仿真程序进行调试、校核和修改；（5）确认仿真程序正确无误后，进行仿真运算；（6）对仿真结果进行分析、评价，并作出决策。25仿真技术的特点：仿真作为一门利用模型进行试验、研究和培训的技术，具有可控、安全、经济、节约时间、允许多次重复的特点。仿真技术的应用领域：仿真技术已广泛地应用于航空航天、航海、国防、原子能、电力、冶金、化工、医学、农业等领域。三、火电厂仿真技术的应用26火电厂计算机仿真技术是一门综合技术。它是以数学、物理、化学、传热学、热力学、流体力学、控制理论、计算机技术、热能动力、电工学、热工仪表及电气仪表等多学科专业的理论为基础，以计算机和各种物理效应设备（它再现真实环境）为工具，对真实火电厂发电机组及其运行进行仿真的技术。27火力发电机组的设备庞大，系统复杂。因此：对设备的研究、对系统的试验、对参数的校正、对运行安全和经济的分析、对运行值班员的培训等，在实际发电机组上直接进行，或是很困难或根本不可能。利用仿真技术的特点、在模型上进行试验研究，在仿真机上培训运行值班员，则上述问题可迎刃而解。

28火电厂仿真技术的应用主要在以下几个方面：对运行值班员进行培训；对设备或系统进行试验、研究、设计、评价、改进、参数优化；对机组运行进行分析、对控制系统的组态和控制规律等进行探讨；对已发故障进行验证，提出反事故措施；对运行规程进行校核。29火电厂应用的仿真机（系统）全范围培训仿真机缩减范围培训仿真机通用型培训仿真机基本原理培训仿真机分析研究仿真机30由于整个计算机仿真过程涉及到：实际系统、数学模型和计算机三个部分。因此，实施有效地仿真必须：

充分认识实际系统、合理建立数学模型、灵巧地应用计算机。这三个方面将是本课程所涉及的内容。31对系统正确的认识是正确的决策的基础，而认识系统必须了解系统的基本特点。系统的特点：（1）系统具有一定功能——不同功能系统，研究方法不完全相同；（2）系统具有相关性——不是元素的简单组合；（3）系统具有整体性——组成部分不可分割；（4）系统是可辨识的——确定的系统才有仿真的可能；（5）系统是有边界的——研究是在一定范围内进行的；（6）系统是可分解的——大系统可分为若干小系统来研究；第三节

系统特点与性质32系统的性质对一个系统进行仿真时，必须明确是一个什么系统、具有什么性质、需作何处理。——对建立系统数学模型极为重要。所谓系统性质是指系统的属性，系统性质通常有：

线性系统非线性系统连续系统离散化系统离散事件系统等之分。33应当明确：（1）不同性质的系统，会有不同的研究方法、不同的仿真方法。（2）描述一个系统时，并非完全按系统的自然特性确定系统性质，而是根据研究的目的、要求、条件等来确定系统模型的类型34第2章

系统数学模型的建立数学模型：——是系统的数学描述，是系统研究的基础，是计算机仿真的依据。352·1建立系统模型的任务（1）确定系统模型的结构

——定义模型性质、确定模型框架和边界条件、明确各环节的特性和相互关系。（2）提供系统模型的数据

——确定系统中各环节特性的定量关系，确定各环节相互间的定量关系（即信号传递的定量关系）。两项任务密切相关，并非各自独立。36注：建模考虑的重点：系统主要的、本质的内容。——数学模型仅仅是系统一种简化﹑抽象的描述，不需要考虑系统的全部细节。实际系统中，不可能只有某种单一的特性（如完全线性、或完全非线性等），但在系统研究时，应将占优势的特性作为重点研究，非优势部分作特殊情况处理；一个系统的数学模型并不是唯一的。但研究目的相同时，研究结果应是一致的；372·2系统模型的分类实物模型：根据相似性建立的形象模型，具有实体性——属物理模拟试验技术范畴。（不介绍）数学模型：用符号和数学方程式表示的系统模型，具有抽象性。（本课程内容）系统模型按描述的状态分按描述的方式分按系统的性质分按求解的方法分按获取的方法分静态模型动态模型连续模型离散模型线性模型非线性模型分析求解模型数字求解模型理论模型黑箱模型反映变量间的相互函数关系反映变量与时间的函数关系是时间的连续函数是时间的离散函数可以用线性微分方程描述不能用线性微分方程描述利于解析法求解利于计算机求解理论推导所得的模型试验研究所得的模型38数字仿真实际上是：——用数值计算技术求解动态数学模型，——系统动态模型在计算机上的数值解。注：（1）一般，仿真所要建立的是系统的动态模型，静态模型包含在动态模型之中。（2）电厂的系统仿真，往往综合采用连续﹑离散﹑线性﹑非线性﹑理论﹑黑箱等模型，以保证模型精度。392·3火电厂模型的特点（1）模型复杂发电厂是一个大系统，它有许多：

1．不同性质的热力设备组成2．不同的能量转换形式3．不同性质与状态的介质4．不同的流动状态5．不同的物理或化学反应6．不同的内扰和外扰7．不同的状态量﹑控制量和操作量8．不同的环节连接方式（串＼并＼反馈＼交叉等）……——设备庞大、系统复杂、参量繁多、交叉影响严重，因而建模工作和描述系统的数学模型也必然十分复杂。数学模型是一个复杂的多元函数组成的数学模型。40（2）线性与非线性并存发电厂有些环节是线性的，有些环节是非线性的。线性特性．汽机功率——调节级压力．抽气量——汽机功率(凝汽式机组)………非线性特性．炉膛传热——温度．流量——压差．转速——油压

………可以应用叠加原理不能应用叠加原理必须有效地协调线性与非线性之间的相互关系必要时,非线性特性——线性化41(3)多为分布参数发电厂生产过程是一个动力循环，工质状态不仅随工况和环境条件变化，而且随时间变化，绝大部分参数是三维空间的函数，具有明显分布参数的特征，

如：过热器的参数(温度)是沿管长、随时间变化的。工质状态（温度）

工况﹑环境﹑条件（过热器管长）

时间建模时多采用分区集总方法，即将三维空间的分布参数简化为一维空间。否则无法求解。42（4）时间常数差别大在发电厂中各设备的动态特性不同，其时间常数（动态响应速度）差别十分悬殊，例如：汽机甩负荷——转速烟温——主汽温时间常数小响应快燃料——汽压减温水量——主汽温时间常数大响应慢建模时应处理好模型的简洁与精确之间的关系对响应速度差异很大的对象，一般分别建立动态数学模型，可采用不同的时间步长进行仿真计算。43（5）模型可分解可集合分解集合把实际系统分成若干子系统，分别建立其模型，然后综合。把若干子系统视为一个整体，统一建立其模型。模型精细﹑复杂模型简洁﹑粗糙因此，在建立复杂系统的动态模型时，应善于对系统进行合理的分解，即按不同的化学、物理过程（如燃烧、传热、流动、做功等），把系统分成许多子系统，这些子系统并不都是简单地按设备的类型来划分的，而根据系统的属性来划分更为合理。在分解系统时，正确处理各子系统间的关系是十分重要的，因此要正确地选择系统中的状态变量，充分保证信息的传递如实际过程那样进行。44（6）模型种类多例如：一台200MW火电机组全范围仿真机，数学方程有6000多个；一台300MW火电机组全范围培训仿真机，数学模型的种类达200多种，构成仿真系统的模型模块达5000多个。为满足实时运行要求，仿真机应采用速度快、容量大的计算机。

（7）模型计算涉及的数据量大数学模型运算涉及的I/O变量、中间变量、常数等数量巨大。例如：一台采用常规仪表盘/台的300MW燃煤发电机组的全范围仿真机有：5000多个I/O变量；10000多个中间变量；

45（8）知识的综合性强这意味着：火电厂建模过程中要综合利用各种知识，解决动力循环中的种种问题。建模是综合知识的体现，它涉及的知识广泛。如：理论知识——物理学﹑化学﹑数学﹑热力学﹑传热学﹑………专业知识——锅炉原理﹑汽轮机原理﹑电机原理﹑机械原理、控制理论﹑………实践知识——结构知识﹑运行知识﹑试验技术﹑………462·4建立数学模型的基本原则本节中，着重介绍建立数学模型的一些共同的、指导性的原则，从总体上了解建模过程中如何处理一些重大问题。而具体的建模步骤和方法，将在后续章中介绍。原则1：必须明确建模的目的和范围：用于研究——用于培训静态模型——动态模型整体特性——局部特性不同的指导思想，所得环节的多少、方程的数目及其变量的多少等都不相同，其差别可能有数十以至成百个之多。47原则2：应建立系统的方框图方框图——用不同的方框来描述系统的各不同部分，各个方框之间依据信号的传递关系连接成一个整体，概括地说明系统的特性。

每个方框——都是由系统的分解而得。

系统的方框图是用来指导系统研究的，它是对系统的最原则的综合。一般来说，可以根据设备的功能、介质的性质和过程扣特点把一个系统划分成许多子系统，子系统又是由许多环节组成的，当不再往下分解时，环节即为分解的极限，从而确定系统的外部边界和内部边界，于是整个框图的雏形便形成了。48系统分解可按：（1）按设备的功能和工作特点分解；如：循环水泵、给水泵、水位控制系统等。（2）按部件的属性分解；如：锅炉各部件属性差别悬殊，可按传热性质相同、热负荷强度相近或介质的性质比较接近来分解，可分解成：炉膛、烟道的烟气侧和工质侧，汽包、下降管和水冷壁，过热器，再热器，减温器，省煤器和空气预热器等。金属的蓄热能力可融合在工质侧之中，不受热部分的联箱可不单独建立环节等。（3）按信号的响应关系分解。在划分系统和环节时，必须明确取得一些什么中间状态响应，最后输出什么响应，从而确定相应环节的方程及其状态变量。系统的输入（如扰动）和输出点，应设置在环节的分界处。当有必要在分布参数的中间提出状态量时，模型应该分成相应的区段，例如同级过热器的高温段和低温段等。系统数学模型建立在每一个方框基础上——分解模型。49原则3：分解模型应具有集合性

建立系统数学模型时，需要进一步考虑能否把一些个别实体组合成更大实体的程度。例：单独研究汽轮机时，为提高精度，可能把它划分成进汽部分、喷嘴、每一级等小实体；而研究发电厂时，汽轮机只是系统中的一个实体，环节就不一定划分太细。

集合是一种简化，在研究系统时，应考虑将某些环节集合处理，从而使模型更为简洁和明确，有利于集中研究某一侧面，突出研究的中心。原则：分解模型不仅要考虑它们集合的可能性，而且能给出它们相互联系的关系式，（该关系式应能反映这个集合的基本特性。）50原则4：信息应具有相关性系统的模型代表一个实际系统的简化形成，它高度概括了系统中变量的相互关系，因此：

模型中应只包括与研究目的有关的信息，模型中应排除与研究目的的无关的信息无关的信息————增加模型的复杂性；增加求解的工作量；容易导致出错。51原则5：应保证模型的准确性这就要求我们根据系统的研究目的和研究范围，把握环节的本质属性，才能做到合理的简化。

简化的程度应由系统研究的精度要求来确定。准确的模型是取得正确仿真结果的基础系统简化与的正确与否密切相关假设条件物理规律依赖于对系统内部化学规律的熟悉程度和概括能力运作机理

522·5电厂仿真建模的一般要求（1）仿真系统所用的数学模型应能够全面描述发电机组在各种可能的过渡过程下的动态性能和所有稳态状态下的稳态性能。因此，发电机组主要系统的动态数学模型应严格地根据质量、能量和动量守恒定律列写基本方程（包括微分方程和代数方程，传热方程、工质热力学状态方程，以及各种状态变量之间的关系式）。是指与热平衡相关的系统。例如主燃料系统、风烟系统、炉膛、主蒸汽系统、汽轮机和抽汽系统、发电机、凝汽器、凝结水系统、给水系统、循环水系统、冷却塔和和加热器疏水系统等。53（2）上述守恒定律可用稳态形式，也可用非稳态形式。只有描述暂态过程的微分方程，在各种状态下其时间常数等于或小于0.5秒时，守恒定律才采用稳态形式。例如：在可压缩流体中的质量平衡方程和低惯性流体（如蒸汽、烟气）中的动量平衡方程。54（3）完整地描述一台发电机组，需要用相关关系式来补充守恒方程。例如：设备特性——诸如泵、阀门的特性可以利用拟合实际数据的相关关系式表达或用表格寻查技术求得。曲线拟合或表格寻查的精度在全部模拟的运行范围内应满足一般精度要求。流体特性——流体特性可以由相关关系式表达或由表格寻查技术求出。流体特性表达式的精度应与模拟的全部运行范围内的稳态精度要求一致。热传递和摩擦损失计算——所采用的相关关系式应当严格地根据公认的工程关系式求出。552·6电厂仿真建模的步骤和方法（1）收集被仿真发电机组的设计资料火电厂仿真机设计的依据是设计资料（不是其物理模型或真实系统）。电站设计资料包括：设备规范、结构设计说明和图纸、性能说明、特性曲线、热力计算书、热力系统管道仪表图（P&ID图）、运行操作规程、报警、保护、操作和控制的逻辑图、控制系统的逻辑图和方框图等。对于已运行的电站，如果对原设计进行了改动，对改动部分应依据改动后的资料。仿真机设计之前尽可能全面地收集到建模所需的资料。56（2）进行初步设计初步设计应利用所收集的主要设计资料、根据对仿真机的要求、按系统和子系统进行。初步设计的主要目的：明确仿真范围，绘制仿真系统图。初步设计内容包括：仿真系统的功能确定及说明、盘/台和DAS（或DCS）系统监视的参数、教练员台监视的参数、系统的简化和假设、故障项目、就地操作仿真项目等。57（3）建立数学模型一次建模（或一次模化）：根据初步设计报告和技术资料，列写系统的数学模型。这样建立的数学模型称为基本模型或初步模型。二次建模（或二次模化）：为了在数字计算机上运算和求解这些模型，必须做进一步简化处理和加工，选择恰当的算法，转换为可用计算机语言编写程序的形式。显然，逐一设备和部件、逐一子系统地列写基本数学模型，必然要花费巨大的工作量。58目前国内外广泛采用“仿真软件开发支持系统”。这些支持系统内都有支持模型软件开发的模块库，库内所有模块都是根据构成火电机组的基本部件和元器件，预先建立好的数学模型，即基本模型（初步模型）,并将它们转换为可在计算机上运行的子程序或算法块，存贮在计算机大容量存贮器的模块库内。在初步设计完成后，在支持软件的支持下，直接利用模块库内的模块和开发工具开发模型软件，并不需要建立基本模型这一过程。592·7热力系统数学模型的表示一、发电厂热力过程的仿真系统，可以认为是一个连续系统。在发电厂中，除了机组的启停和甩负荷等特殊工况外，整个生产过程是一个能量转换的连续过程，因此，反映发电厂热力过程的仿真系统，可以认为是一个连续系统。二、连续系统动态特性的数学模型通常用高阶微分方程来描述。由于方程的可转换，反映连续系统动态特性的数学模型又可以具体地表示为：高阶微分方程、一阶微分方程组、传递函数、复变函数、状态方程。60三、发电厂热力系统的仿真，既可采用连续系统的数学模型，也可采用离散系统的数学模型。既可采用线性系统的数学模型，也可采用非线性系统的数学模型。由于发电厂热力系统是一个大系统，其设备组成非常复杂，有线性模型，也有非线性模型。如果采用非线性模型进行系统仿真，可以直接根据系统的实际情况建立非线性模型。只要能建立系统的非线性模型，就有可能简化成线性模型进行系统仿真。采用线性模型的优点是通过线性代数的运算，能很快地得到所需要的解，但由于建立非线性模型时已对系统作了不少假设或简化，在线性化过程中，还要作进一步简化，所以模型的精度并不是总能保证的。61四、热力系统仿真广泛采用了黑箱模型。原则上，理论建模法是比较严谨的，但由于系统复杂，某些机理一时弄不清楚，或模型过于复杂以致得不到所需要的简化形式，使我们不得不借助于实验数据，用有限的变量来建立系统的所谓黑箱模型。正是由于这些原因，热力系统仿真中，广泛采用了黑箱模型。发电厂热力系统仿真中的建模方法与控制系统仿真中的建模方法是不完全相同的，它往往是线性模型、非线性模型、理论模型、黑箱模型等的综合采用，——以保证模型的精度。622·8连续系统的数学模型表达方式

（1）以高阶微分方程表示的规范式一个连续系统可用高阶微分方程表示如下：自由项强迫项(2-1)式中：y=y（t）——系统的输出量

u=u（t）——系统的输入量

ai，cI———系数63当引入算子时，则式（2-1）可表示为：pny+a1pn-1y++an-1py+any=c0pn-1u+c1pn-2u++cn-1u或：（其中a0=1）(2-2)64（2）以一阶微分方程组表示的规范式一个连续系统一般是由一些局部环节组成的，有时我们不把它们综合成系统的高阶微分方程，而直接以一阶微分方程组来求解，也是很方便的。其表达式如下：aij随时间而变——时变系统aij不随时间变——定常系统Ui(t)全部恒等于零——齐次线性微分方程组Ui(t)有一个不为零——非齐次线性微分方程组式中：yi——局部输出量；ui——局部输入量；aij,——系数（2-3）65（3）高阶微分方程与一阶微分方程组的关系高阶微分方程与一阶微分方程组可互为转化。例如，已知一个系统的n阶微分方程为：令

y1y2y3…

yn（2-4）式中：y1，y2，，yn为自变量t的n个未知函数。求解式（2-4）相当于求解式（2-5）(2-5)66（4）传递函数表示的规范式如果式（2-1）高阶微分方程y和u的各阶导数（包括零阶）的初值均为零——各初始条件全为零时，对式（2-1）的两边进行拉普拉斯变换有：（2-6）即：sn—替换—dn/dtns为拉普拉斯算子

Y(s)—替换—y(t)Y(s)为输出量y(t)的拉普拉斯变换

U(s)—替换—u(t)U(s)为输入量u(t)的拉普拉斯变换

此式与式（2-2）比较可知，在初值为零的情况下，用算子P表示的式子与用传递函数G(s)表示的式子在形式完全相同。（2-7）将上式整理可得系统传递函数G（s）的规范式为：：67（5）复变函数表示的规范式一个系统也可用复变函数描述。如果用jω替换传递函数中的s，对应于式（2-7）的系统，其复变函数表达的规范式为：（2-8）68（6）状态方程表示的规范式A、问题的提出：上述模型只描述了系统输入与输出之间的关系未描述系统内部的情况称为：外部模型仿真在计算机上对系统的数学模型进行试验必须在计算机上复现（实现）系统即复现系统输入量复现系统输出量复现系统的内部变量又称：状态变量69仿真要求将系统辨识或其它方法建立的外部模型转化为：内部模型——状态空间模型——状态方程控制理论中称之为：实现问题70B、状态方程的基本概念设传递函数表示的多变量系统为：传递函数G(s)U1（s）U2（s）UL（s）Y1（s）Y2（s）Ym（s）当引入状态变量后，框图变成了：状态变量xi

UYlnm71①

L个输入变量uI，输入变量的集合可用输入向量U来表示；系统中共有三种变量，即：——U、Y、X均为时间的函数。（2-9）（2-10）②m个输出变量yI，输出变量的集合可用输出向量Y来表示；（2-11）③n个状态变量xI，状态变量的集合可用状态向量X来表示：72且在t时刻，有：状态向量可表示为：输出向量可表示为：（2-12）（2-13）系统的状态方程系统的输出方程式中：X（t0）——初始状态向量

U（t0,t）——初始时刻t0到t的系统输入向量

f[]——表示一个单值函数

g[]——也表示一个单值函数若用微分方程表示，则其表示式为：（2-14）（2-15）若用线性微分方程表示，则其表示式如下：73（2-16）（2-17）系统数学模型的状态方程规范式在状态空间中：A（t）——(nn阶)系统矩阵（系数矩阵、状态矩阵）B（t）——(nl阶)控制矩阵（驱动矩阵）C（t）——(mn阶)输出矩阵D（t）——(ml阶)传递矩阵（传输矩阵）当A(t)、B(t)、C(t)、D(t)随时间而变化时，系统称为变系数系统（时变系统）；当A(t)、B(t)、C(t)、D(t)不随时间而变时，可用A、B、C、D来表示，系统则称为常系数系统（定常系统）。D(t)U(t)项表示系统的输入和输出通过D(t)直接联系的部分。在普通的控制系统中，通常D(t)=0，则式（2-17）变成：（2-18）74C、状态方程规范式说明1、某一时刻t时，系统的内部状态Ｘ(t)可用n维空间中的一个点来描述，随着时间的变化，X(t)在这个空间中描绘出连续的轨迹。————因此，在任意时刻，系统的内部状态是由状态空间中X(t)的位置决定的，这个位置可按状态方程的解求得。2、输出向量Y(t)，可用矩阵和向量的简单运算不难求得。75D、系统状态方程的建立当高阶微分方程转换为状态方程时，视强迫运动项（等式的右边项）是否具有导数项而有很大的差别。下面以单输入/单输出系统为例，介绍几种建立状态方程的基本方法：a、微分方程的强迫项无导数项时的状态方程（规范式Ⅰ）在该情况下，式（2-1）中的微分方程为：(2-19）其中：今引进n个状态变量一阶微分方程组76把上述一阶微分方程组写成矩阵形式，可得：（2-20）(2-21)77令：则可得到以状态方程表示的规范式（I）（2—22）（2—23）显然，状态变量的初值为：78注：状态方程的规范式所包含的内容为与所描述系统相对应的：状态方程X输出方程Y（输出向量Y）状态变量xI（或状态向量X）状态矩阵A控制矩阵B输出矩阵C传递矩阵D输入向量U状态变量初值xI（0）等——（下同）.79[例]：已知描述某系统的微分方程为：且：试求该系统的状态方程，并将给出的初始条件用状态变量的初值表示。[解]：令则有将上式写成矩阵形式，80相应地输出方程为：状态变量的初值为：即可得状态方程：规范式（Ⅰ）81此时：b、微分方程强迫项有导数项时的状态方程（规范式Ⅱ、Ⅲ、Ⅳ）当引入算子时，则上式可表示为：pny+a1pn-1y++an-1py+any=c0pn-1u+c1pn-2u++cn-1u或：（其中a0=1）82Ⅰ、当初始条件按引入的各状态变量的初值给出时：设上式中的u为：引进n个状态变量：即：pjx=xj+1

(j=0,1,2,,n-1)（2-24）则式（2-24）可写为：83或：取a0=1,则有：（*）将式（2-24）代入式（2-2）可得：（**）由式（*）和（**）可得：（2-25）——规范式（Ⅱ）84状态变量的初值为：其中：规范式（Ⅱ）的应用：是系统的数学模型可直接写成状态方程形式，且初始条件也是按各状态变量的初值给出时，计算将十分方便。85Ⅱ、当初始条件按输入u和输出y及其各阶导数给出时此时，若将微分方程转变为状态方程描述系统，则必须将给出的初始条件转换成各状态量的初值。假定给出的微分方程组为：（2-30）即u和y的导数阶次相等若令：则有：若将式（2-30）两边同时取不定积分一次并令：则有：86若将式（2-30）两边同时取不定积分二次并令：则有：…………同理有或：则有：即：（j=0，1，2，……n）（2-31）（#）又因为在j=0时，由式（#）有：则在a0=1时，有：求初值用87此时由式（#）可得：由上述讨论，可得式（2-30）在a0=1时的状态方程和输出方程如下：88（2-32）（2-33）——规范式（Ⅲ）89当c0=0（即：输入u的导数阶次＜输出y的导数阶次）时，规范式（Ⅲ）可表示为：（2-32’）（2-33’）——规范式（Ⅲ’）90状态变量的初值：当高阶微分方程给出输入量u和输出量y及其各阶导数的初始值时，代入式（2-31），即可求出各个状态变量的初始值。即：其中：(j=1，2，……n)91[例]：已知描述某系统的微分方程为：式中：u为单位阶跃函数，且初始条件为：试求该系统的状态方程和输出方程，并将给出的初始条件用状态变量的初值表示。[解]：根据题意，有：n=3，a0=1，a1=2，a2=3，a3=4，c0=0，c1=0，c2=5，c3=6可由规范式（Ⅲ’），得系统的状态方程和输出方程：92其中，各状态变量由式（2-31）可求得：93状态变量的初值由式：可求得：应该指出：一个系统可以引入不同组合的状态变量。那么，所设的状态变量不同，获得的状态方程也会不同。即一个外部模型，可能有很多不同的内部模型（状态方程）——一个系统的状态方程形式并非是唯一的。94Ⅲ、当初始条件按输入u和输出y及其各阶导数给出时，

且具有唯一解的状态方程——规范式（Ⅳ）：（j=0，1，2，……n）（2-31）由式（2-31）可知，所设的各状态变量与u的导数项有关，因此而得的状态方程:（2-32）（2-33）——规范式（Ⅲ）其解可能不是唯一的。95设系统的数学描述为：要得到唯一解，必须设置一组状态变量，不包括u的导数项。现代控制论为此提供了另一种变换法。（2-34）如取以下n个变量作为一组状态变量：（3-35）且式中，0，1，2，，n分别为：a0=196（2-36）则可保证状态方程解的唯一性。根据上面一组状态变量，对于式（2-34）所描述的系统，可得出如下规范式：97（2-37）（2-38）式中：状态方程规范式（Ⅳ）98几点结论：1．对于以上各规范式：若A，B，C，D随时间变化——称为：变系数系统或时变系统若A，B，C，D不随时间变——称为：常系数系统或定常系统2．一个系统的数学模型并非唯一的。表示一个系统的模型，可用不同的数学方法。最常用的方法为传递函数，其次是微分方程，状态方程。3．采用微分方程进行连续系统的仿真时：微分方程一阶微分方程组代数方程采用变量分离法一般把转换成转换成994、采用状态方程进行连续系统的仿真时：（a）应根据初始值是否为零，选择相应的规范式。初始值为零时可选用规范式（Ⅰ）或规范式（Ⅱ）。（b）若高阶微分方程的初始条件是输入u，输出y及其导数的初值，应变换成状态方程的初值。（c）若微分方程的强迫项含有导数项，可采用规范式（Ⅱ）或规范式（Ⅲ），但此时可能没有唯一的解，若要获得唯一的解，则应采用规范式（Ⅳ）。5、对已确定的系统，仿真所采用的数学方法可以不同，但研究结果应与系统的实际情况一致。100

2.9系统非线性模型的线性化

在很多时候我们都用线性模型，因为实际当中的非线性模型很难求解，因此，我们用“等效”的线性方程来代替。线性系统的优点：求解容易、具有可叠加性。非线性线性化的条件：变量对预定的工作点偏离很小，即认为系统的变量在动态过程中偏离工作点不远。101非线性方程y=f(x)在（x0,y0)处展开得：

y=f(x)=f(x0)+df(x-x0)/dx+d2f(x-x0)2/2!dx2+······由于x-x0很小，所以(x-x0)2，(x-x0)3可忽略。于是有：

y=f(x0)+df(x-x0)/dx

y=y0+df(x-x0)/dx

△y=K△X102特点：1、在工作点小范围变化时，不用此模型，不同的工作点线性模型不同。2、线性程度越好的模型其允许的波动范围愈大。3、有些模型存在折点、断点等情况时，不能进行这样的处理。4、进行线性处理会引入误差。103第3章面向微分方程的数值积分法仿真104数值积分是数值分析的一个基本问题。也是复杂计算问题中的一个基本组成部分。数值积分往往用极简单的方法就能较好地得出对所求解的具体数值问题的解答。但数值积分的难点在于计算时间有时会过长，有时会出现数值不稳定现象。另外，数值积分的理论性较强。其理论和方法都已经比较成熟，计算精度也比较高。1053.1仿真中研究数值积分法的意义数值解的一种近似方法。对于连续系统的高阶微分方程，可化为若干个一阶微分方程组成的方程组。数值积分法是求解微分方程：例如：下式所示的状态方程可以化为一个一阶微分方程组（3-1）（3-2）106所以，连续系统的仿真就是从给定的初始条件出发，对描述系统动态特性的常微分方程或常微分方程组进行求解，从而得到系统在一定输入作用下的变化过程。在求解这些微分方程时，最常用、也是最有效的一种方法就是数值积分法。1073.2数值积分法仿真的基本原理对微分方程（3-1）两端同时取积分，可得当时，上式变为：（3-4）（3-5）将积分项拆成两项（3-6）108故上式可写为：（3-7）——此式是方程（3-1）在tn+1时刻的精确解。

在数值解法中，希望用近似解：代替准确解，其中：（3-8）为为为的近似值令：——称为计算步长或步距式（3-8）是从常微分方程（3-1）出发建立的离散数学模型——差分方程。为为为为为为109由此可见，数值积分法就是在已知微分方程初值的情况下，求解该方程在一系列离散点处的近似值，其特点是步进式——根据初始值逐步递推地计算出以后各时刻的值。从式（3-8）可知，数值积分法的主要问题归结为如何对f（t，y）进行数值积分——求出f（t，y）在区间[tn，tn+1]上定积分的近似值Qn。采用不同的方法求Qn，就出现了各种各样的数值积分方法。不同的数值积分将对求解的精度、速度和数值稳定性会产生不同的影响，这将在下述内容中具体介绍。110数值积分法种类繁多，在此从实用角度介绍几种基本的方法3.3欧拉(Euler)法3.3.1简单欧拉法

欧拉法是一种最简单的数值积分法，对于方程：在区间[tn，tn+1]上求积分，得到：若区间[tn，tn+1]足够小，则[tn，tn+1]上的f（t，y）可近似地看成常数f（tn，yn）。故可用矩形面积近似代替111即：tntn+1f（tn，yn）于是有：（3-9）将此式写成差分方程为：（3-10）——著名的欧拉公式1123.3.2改进的欧拉法

如果用梯形面积而不是矩形面积来代替每一个小区间上的曲线积分，就可以提高计算精度，梯形法的计算公式为：（3-11）式中的右端含有待求量yn+1，因而它是隐函数形式。这种方法不能自行启动运算，需要依赖其它算法的帮助。每次计算都用欧拉法算出y（tn+1）的近似值，以此计算近似值，然后利用梯形公式（3-11）求出修正后的。即有：帮助方法：113预估式校正式（3-12）——改进的欧拉公式3.3.3几个基本概念

简单的欧拉法是用前一时刻tn的yn求出后一时刻的yn+1，这种方法称为单步法，它是一种自行启动的算法。如果求yn+1时需要tn,

tn-1,

tn-2……时刻yn,

yn-1,

yn-2……的值，这种方法为多步法（改进的欧拉法为两步法），它是一种不能自行启动的算法。1、单步法与多步法114简单的欧拉法表达式的右端，计算所用的数据均已求出，这种公式称为显式公式。改进的欧拉法表达式的右端，有待求量，这种公式称为隐式公式。隐式公式不能自行启动，需要用预估-校正法。单步法和显式在计算上比多步法和隐式方便，但有时为了满足精度、稳定性等方面的要求，需要采用隐式算法。2、显式与隐式1153、截断误差这里用泰勒级数为工具来分析数值积分公式的精度假定yn是精确的，用泰勒级数表示处的精确解，即：显然，简单的欧拉法是从以上精确解中取前两项之和来近似计算，每一步由这种方法引入的误差称为局部截断误差，简称截断误差。简单的欧拉法的截断误差为：不同的数值方法有不同的截断误差。一般若截断误差为，则方法为r阶的。所以方法的阶数可以作为衡量方法精确度的一个重要标志。116由此可见，简单的欧拉法是1阶精度；改进的欧拉法由于采用了平均斜率，相当于取了泰勒级数的前3项，因此为2阶精度。

分析欧拉法截断误差的思想，同样也适用于其它数值积分方法。4、舍入误差由于计算机的字长有限，数字不能表示得完全精确，在计算过程中不可避免地会产生舍入误差。舍入误差与计算步长成反比。如果计算步长小，计算次数多，则舍入误差就大。舍入误差除了与计算机字长有关以外，还与计算机所使用的数字系统、数的运算次序以及计算所用的子程序的精度等因素有关。1175、数值解的稳定性问题采用数值积分法求解稳定的常微分方程，应该保持原系统的稳定特征。但是：（1）在计算机逐步计算时，初始数据的误差及计算过程的舍入误差对后面的计算结果将产生影响。（2）如果计算步长取的不合格，有可能使仿真出现不稳定的结果。下面我们简单讨论一下这个问题。差分方程的解与微分方程的解类似，可分为特解和通解两部分。与稳定性有关的是方程的通解，它取决于差分方程的特征根是否满足稳定性条件。例如，为了考查欧拉法的稳定性，我们研究检验方程（TestEquation）：其中，为方程的特征根118显然，要使该差分方程是稳定的，必须使下式成立。即：表明：为使数字仿真稳定，对计算步长应有所限制。另外，稳定性还与系统的特性以及数值积分方法有关。上述分析欧拉法稳定性的思想，同样适用于其它数值积分方法。（3-14）（3-15）1193.4龙格-库塔（Runge-Kutta）法由前面的分析可知，将泰勒展开式多取几项以后截断，就能提高截断误差的阶数和计算精度。然而，直接采用泰勒展开方法要计算函数的高阶导数，运用起来不便。龙格-库塔方法的基本思想是：用几个点上的函数值的线性组合代替函数的各阶导数，然后按泰勒级数展开确定其中的系数，这样既可避免计算高阶导数，又可提高积分的精度。龙格-库塔法有多种形式，以下从实用角度直接给出公式的形式和相应的精度。3.4.1龙格-库塔方法的基本思想1203.4.2二阶龙格-库塔方法2阶龙格-库塔方法的公式为：（3-16）上式表示的数值解是用的泰勒级数在2阶导数以后截断所求得的，因此称为2阶方法。故2阶龙格-库塔法与式改进的欧拉法相比，实质完全相同。所以改进的欧拉法实质上是2阶龙格-库塔法。截断误差为：121实时仿真的2阶龙格-库塔方法（3-17）其截断误差为：1223.4.3四阶龙格-库塔方法4阶龙格-库塔法是一种最常用的方法。其经典表达式为：（3-18）其截断误差为：123显然，4阶龙格-库塔法的计算量较大，但计算精度较高，在比较不同算法的计算精度时，常以它的计算结果作为标准。实时仿真的4阶龙格-库塔方法（3-19）124以上公式都是标量形式，如果要换成向量形式，只要把式中的标量y，f，k换成向量Y，F，K即可。从理论上讲，可以构造任意阶数的龙格-库塔方法，但是，精度的阶数与计算函数值f

的次数之间并非等量增加的关系，见下表所列：

由此可见，4阶经典龙格-库塔方法有其优越性，而4阶以上的龙格-库塔方法计算f所需的次数比阶数多，增加了计算量，从而限制了更高龙格-库塔方法的应用。对于大量的实际问题，4阶方法已可满足精度要求，所以得到了广泛的应用。125我们仍采用检验方程进行讨论，对它利用泰勒级数展开得：3.4.4龙格-库塔公式的稳定区域对于，有，将它代入上式得：（3-20）（3-21）令：——将其代入上式得到该式的稳定条件为：126（3-22）由此稳定条件，下表给出了各阶龙格-库塔公式的稳定区域。表3.2龙格-库塔公式的稳定区域在使用龙格-库塔公式时，选取的步长应使落在稳定区域内。否则，在计算时会产生很大的误差，从而得不到稳定的数值解。127例如：用4阶龙格-库塔方法解：（1）用解析法求解：本例是稳定的；（2）用数值法求解：当h=0.1时，数值解是稳定的;当h=0.2时，数值解就不稳定了，这时因为：此数值在稳定区间以外，所以数值解不能收敛。这种对步长有限制的数值积分法称为条件稳定积分法。128从4阶龙格-库塔方法的稳定条件：中可以看出，系统的特征根越大，需要的积分步长就越小，这一点可作为选择步长的依据。步长的大小，除了与数值积分方法的阶数有关外，还与方程本身的性质有关。除以上介绍的欧拉法、龙格-库塔法外，数值积分方法还有许多种，如亚当姆斯方法、吉尔方法等等。此处不一一介绍。1293.5计算方法的选择数值积分方法很多，在实际使用时存在一个选择问题。对于一个具体问题，如何选择具有一定的难度，至今尚无一种确定的方法。一般来说，数值积分方法的选择应考虑的因素有：1、精度要求数值积分法的精度受截断误差舍入误差积累误差的影响，而这些误差与积分方法、步长、计算时间、计算机精度等有关。130一般：（1）积分方法的阶数越高，截断误差越小，精度越高；（2）步长越小，精度越高；（3）多步法的精度高于单步法；（4）隐式算法的精度高于显式算法。因此，当需要高精度时，可采用高阶、多步、较小步长、隐式算法。但是，步长的减小往往会增加迭代次数并增大舍入误差和积累误差。与此相反，在精度要求不高时，最好使用低阶方法。131计算速度取决于计算步数、每步积分所需的时间。而每步的计算时间又与积分方法有关，它主要取决于计算导数的次数（4阶龙格-库塔方法每步要计算4次导数），导数的计算是最费时的；为了加快计算速度，在积分方法已定的条件下，应要保证精度的前提下，尽量选择较大步长，以缩短积分时间。2、计算速度132数值解的稳定性必须保证，否则，计算结果将失去真实意义。从稳定性来看，不同的数值积分算法有不同的稳定性。应用时应控制步长h，使数值积分算法在稳定域内。3、数值稳定性133步长的选择很重要，步长过大会增大截断误差，甚至出现数值不稳定现象，过小了又因增加了步数，而使舍入误差增大。所以，仿真的总误差与步长的关系不是单调函数关系，而是一个具有极值的函数。如图所示。3.6计算步长的选择误差总误差截断误差舍入误差步长由图可知，存在一个最佳步长。134当积分方法确定以后，在选取步长时，需要考虑的一个重要因素就是仿真系统的动态响应特性。如果系统的动态响应快，导数变化激烈，则应取高阶方法和小步长进行计算。为了保证数值稳定性，步长应限制在最小时间常数（相当于最大特征值的倒数）的数量级上。用经验方法选取步长的两种方法：（1）由系统方程中最小时间常数Tmin来决定：（2）由系统开环频率特性的剪切频率来决定：1353.7面向微分方程的仿真程序构成一般的仿真程序的组成可用下图表示：输入或设置参数子程序主程序运行管理子程序解答存储及打印子程序在线绘图子程序仿真系统计算子程序图中各方块的功能为：主程序：实现对整个仿真计算的逻辑控制。输入或预置参数模块：输入系统参数初值，计算步长、时间等参数。运行管理模块：这是数字仿真程序的核心，对仿真计算进行时间控制，以保证计算机按要求进行计算及输出。计算模块：根据被仿真的系统及所选的仿真方法编写的计算程序。输出及显示模块：将仿真结果以数据或图表的形式输出给用户。1364阶龙格-库塔法仿真程序设计

目前，应用于系统仿真的商用软件包诸多，但无论是开发者还是应用者，了解程序的设计思想、分析程序的结构与功能，完善和编写仿真程序都是必要的。因此，此处通过对经典龙格-库塔方法（即式（7-22））仿真程序的介绍，使读者了解和掌握编写仿真程序的基本技术。为了便于说明程序的实质，我们用具体的例子来讨论。假设有一个2阶系统：令：则：137以下是用C语言编写的原理性程序：

程序中变量说明：i，n: 分别为循环控制和积分器个数变量。h，t，tmax：分别为步长、时间和仿真总时间变量。yn[10]: 存放tn时刻y值的临时变量数组。y[10]: 存放y值的变量数组。doty[10]: 存放的数组。k1[10]，k2[10]，k3[10]，k4[10]:分别为存放K1，K2，K3，K4的数组。程序清单：（“/*”和“*

/”之间的文字是非执行语句，作为注释使用。）#include<stdio.h>voiddiffeq();inti,n;floath,tmax,t,yn[10],y[10],doty[10],k1[10],k2[10],k3[10],k4[10];main() /*主函数*/138{n=2; /*设置积分器个数*/y[1]=0;y[2]=0 /*设置积分器初值*/h=0.01;tmax=20;t=0/*设置步长、仿真总时间和初始时间*/whilt(t<=tmax) /*判断是否仿真时间到*/{printf(“t=%f”,t)； /*打印时间*/printf(“y[1]=%fy[2]=%f\r\n”,y[1],y[2]);/*打印积分器输出*/for(i=1,i<=n,i++) yn[i]=y[i]; /*保留tn时刻的初值*/dif(); /*计算*/for(i=1;i<=n;i++){k1[i]=doty[i]; /*计算K1*/y[i]=yn[i]+(h/2)*k1[i];}139t=t+h/2dif(); /*计算*/for(i=1;i<=n;i++){k2[i]=doty[i]; /*计算K2*/y[i]=yn[i]+(h/2)*k2[i];}dif(); /*计算*/for(i=1;i<=n;i++){k3[i]=doty[i]; /*计算K3*/y[i]=yn[i]+(h/2)*k3[i];}t=t+h/2dif(); /*计算*/140for(i=1;i<=n;i++){k4[i]=doty[i]; /*计算K4*/y[i]=yn[i]+(h/6.0)*(k1[i]+2*k2[i]+2*k3[i]+k4[i]); /*计算tn+1时刻的y值*/}}}voiddif() /*计算变量导数的函数*/{floatu;u=1; /*计算u=1(t)*/doty[1]=y[2];doty[2]=4*u-2.828*y[2]-4*y[1];}下面是程序的主函数main()的处理流程图

141四阶龙格-库塔方法的程序处理流图142第4章

连续系统仿真的离散相似法“离散相似法”——将一个连续系统进行离散化处理,然后求得与它等价的离散模型(差分方程)的方法.143获取离散相似模型的两个途径:

（1）对传递函数作离散化处理得离散传递函数——称为“频域离散相似模型”；（2）基于状态方程离散化——称为“时域离散相似模型”；4.2状态方程的离散化4.1传递函数的离散化本章介绍:4.3典型环节的离散化模型4.4常用非线性环节的仿真4.5离散相似法仿真的特点1444.1传递函数的离散化对连续系统进行数字仿真可以先在系统加入虚拟的采样器和保持器，如图4-1所示，然后利用Z变换的方法求出系统的脉冲传递函数，再从脉冲传递函数求出对应于系统G(s)的差分方程。保持器TTuy图4-1连续系统离散化结构图

145如此处理的理由是：（1）由采样定理可知：当采样频率和信号最大频率满足：的条件时，可由采样后的信号唯一地确定原始信号。

（2）把采样后的离散信号通过一个低通滤波器，即可实现信号的重构。常用的滤波器是零阶保持器和一阶保持器等。

因此，研究系统G(s)，在一定的条件下可能等效地研究图4-1所示的系统。注意：图4-1所示系统的采样开关和保持器实际上是不存在的，而是为了将(4-1)式离散化而虚构的。146根据图8-1,有脉冲传递函数:其中Gh(s)是保持器的传递函数。若选择不同的保持器，则可得不同的G(z)，见表4-1。（4-1）保持器的传递函数Gh(s)

脉冲传递函数G(z)零阶：

一阶：

三角形：

表4-1不同保持器的G(z)147表8-1中的G(z)还要求对或等进行Z变换，这样做有时不太方便，如果能对G(s)直接进行Z变换就方便多了。

因此，对表8-1中的各式还可以进一步加以简化，得出二次加入虚拟采样器和保持器的脉冲传递函数。例如:

当采用零阶保持器时需要对进行Z变换，它相当于G(s)与1/s串联，如图4-2所示的模型。TT图4-2若在积分环节之前再加一个采样器和保持器，如下图:TT图4-3T148那么就可以对G(s)及1/s分别求出脉冲传递函数，然后相乘，即:

因此,二次加入采样器、零阶保持器后的脉冲传递函数为:（4-2）（4-3）

从上述可看出，虚拟采样器和保持器可以不止一次地使用。在必要的地方加入，可为求脉冲传递函数带来方便。但须注意，由于保持器的频谱特性并非理想矩形，所以每加一次采样器和保持器都会带来误差。因此(4-3)式较之表4-1中的式子误差要大些，所以要尽量减少虚拟采样器和保持器的使用。

149下面举例说明其使用方法。若:则根据表4-1，当加零阶保持器时，可得:所以得差分方程为:也可根据(4-3)式来求G(z)：

其对应差分方程为:（4-4）（4-5）150以上(4-4)式或(4-5)式是按一次和二次加入虚拟采样器和保持器时分别求得的差分方程。后者较前者误差要大些。有了上面的知识，不难将其推广应用到结构图表示的系统，求其等效的差分方程。在此不再多叙。1514.2状态方程的离散化假设连续系统的状态方程为:

若人为地在系统的输入端及输出端加上采样开关，同时为了使输入信号复原为原来的信号，在输入端还要加一个保持器，如图所示。

保持器TTux假定为零阶保持器:——输入向量的所有分量在任意两个依次相连的采样瞬时为常值。例如:对第n个采样周期u(t)=u(nT)。

T为采样间隔(周期)（4-6）图4-4采样控制系统结构图

152若对方程(4-6)式两边进行拉普拉斯变换，得：故对(4-7)式反变换可得：即：以(sI-A)-1左乘上式的两边可