matlab对线性系统稳定性的分析

MATLAB对线性系统稳定性的分析

摘要：本文对线性系统从时域、复域和频域进行了稳定性分析，总结了控制系统的主要判据，分析过程简单，结合实例验证了其真实性、有效性。关键词：线性系统稳定性MATLAB

引言：一个控制系统要能正常工作,必须首先是一个稳定的系统,即当系统受到外界干扰后,虽然它的平衡状态被破坏,但在外扰去掉以后,它仍有能力自动地在平衡状态下继续工作。在已知一个系统的系统函数或状态空间表达式时,就可以对其系统的稳定性进行分析。但当系统的阶次较高时,绘图和计算需要花费大量的时间和精力。MATLAB是一套高性能的数值计算和可视化软件,并拥有几十个工具箱,借助MATLAB的系统工具箱,就可以直观、方便地分析系统的稳定性。

1、控制系统稳定性定义

关于稳定性的定义有许多种，较典型的说法有两种：一种是由俄国学者李雅普诺夫首先提出的平衡状态稳定性，另一种指系统的运动稳定性。对于线线控制系统而言，这两种说法是等价的。根据李雅普诺夫稳定性理论，线性控制系统的稳定性可以定义如下：若线性控制系统在初始扰动的影响下，其过渡过程随着时间的推移逐渐衰减并趋向于零，则称该系统为渐近稳定，简称为稳定；反之，若在初始扰动影响下，系统的过渡过程随时间的推移而发散，则称系统为不稳定。由上述稳定性定义可以推知，线性系统稳定的充分必要条件是：闭环系统特征方程的根都具有负实部，或者说闭环传递函数的极点均位于左半S开平面（不包括虚轴）。

2、系统稳定性分析方法概述

在经典控制理论中，常用时域分析法、复域分析法或频率分析法来分析控制系统的性能。不同的方法有不同的适用范围，下面对上述方法进行具体研究。

2.1时域分析法

在经典控制理论中，时域分析法是一种直接在时间域中对系统进行稳定性分析的方法，具有直观、准确的优点，并且可以提供系统时间响应的全部信息。在时域分析系统的稳定性，必须研究在输入信号作用下，当时间t趋于无穷时，系统的输出响应趋于最终期值h（∞）。显然，一个稳定的系统，其时域响应曲线必须是衰减的。

2.2复域分析法

在复域中进行系统稳定性分析，尤其当系统参数K 的变化时，选定合适的参数范围使系统达到所需要稳定要求。有两种方法：一是直接法，即对于较易得到系统闭环传递函数的场合，直接求出系统所有闭环极点，判断是否都具有负实部来确定系统的稳定性；二是根轨迹法，利用系统开闭环传递绘制根轨迹，由线性系统稳定的充分必要条件：闭环传递函数的极点均位于左半S 开平面（不包括虚轴），确定使根轨迹在左半S 开平面部分时参数范围为系统稳定的区域。

2.2.1直接法

假设闭环传递函数为)(s φ=)

()(s D s M ，则其特征方程0)(=s D 写成一般形式： 0...)(122110=+++++=---n n n n n a s a s a s a s a s D ()00>a (1)

若n ≤2，可直接求取其特征方程根（即闭环极点）来判断系统稳定性，即使（1）有待定参数，也容易求出特征方程根的一般形式，但对于求取n>3的高阶系统特征方程式的根很麻烦，所以对高阶系统一般都采用间接法来判断稳定性，在时域中常采用间接方法是代数判据（也称劳斯判据）。

2.2.2根轨迹法

根轨迹法是一种图解方法，这种方法是根据系统开环零、极点的分布来研究系统中可变参数变化时，系统闭环特征根的变化规律，从而研究系统的稳定性。 因此，根轨迹法在控制系统的分析和设计中是一种很实用的工程方法。它的最大特点是能够很清晰地了解到闭环特征根的分布，一目了然地得出系统稳定时参数的取值范围，并且不必求出系统的闭环传递函数，适用于较复杂系统。根轨迹法的关键环节就是能够正确地绘制出系统的根轨迹，简单根轨迹可用试探法绘制，复杂根轨迹则应利用其绘制基本规则进行绘制。

2.2.3频域分析法

频域分析法是应用频率特性研究系统的一种经典方法，以系统的频率特性为数学模型， 用bode 图或其他图表作为分析工具。 当系统的开环传递函数表达式不易求出，就无法应用代数判据或根轨迹法判断闭环系统的稳定性，此时应用频率稳定判据就非常方便。其前提条件就是要正确地把系统的频率特性绘制成曲线，常用的频率特性曲线大致有三种：幅相曲线（极坐标图）；bode 图，也称为对数频率特性曲线；对数幅相曲线（尼科尔斯图）。曲线的绘制可根据系统的开环频率特性的表达式通过取值描点法、叠加法绘制根轨迹草图，或利用MATLAB 等计算机辅助工具来实现。

3、MATLAB 实现系统稳定性分析

3.1时域分析法判断系统的稳定性

系统模型为WK (S)=)2)(550

++s s （，单位负反馈。利用MATLAB 工具箱提供的时域响应函数， 给该系统施加单位冲激，观察它的响应，分析稳定性。

程序如下：

程序中num 为开环传递函数分子系数矩阵，den 为分母系数矩阵。系统的稳定性是指系统在遭受外界扰动偏离原来的平衡状态，当扰动消失后，系统自身仍有能力恢复到原来平衡状态的一种能力。从图1可以很直观地看出该系统是稳定的。

图1 单位冲激响应图

3.2直接判定法

根据稳定的充分必要条件判别线性系统的稳定性， 最简单的方法是求出系统所有极点，并观察是否含有实部大于0的极点，如果有，系统则不稳定。比如, 要判断系统函数2

22)(2+++=s s s s H 的系统是否稳定,可利用MATLAB 快速求出其零极点并绘出零极点的分布图,程序如下:

MATLAB 运行结果如下:

从运行结果看,系统极点P 值实部全部为负,得知极点全部位于S 左半平面,可判断该系统为稳定系统。

对于简单的系统函数可以通过数学运算直接求得极点来分析系统稳定性,然而实际的控制系统大部分都是高阶系统，这样就面临求解高次方程，求根工作量大但在MATLAB 中只需调用den 函数即可，这样就可以由得出的极点位置直接判

定系统的稳定性。比如系统函数为1

22532423)(2345234+++++++++=s s s s s s s s s s H ，程序如下:

MATLAB 运行结果如下:

图2 系统零极点分布图(×表示极点, o 表示零点。)

从计算机结果以及零极点图2可以看出,该系统的极点并不都在s 左半开平面,有一对共轭极点位于S 右半开平面,所以该系统不稳定。

3.3轨迹法判断系统的稳定性

MATLAB 控制工具箱中提供了rlocus 函数， 来绘制系统的根轨迹，利用rlocfind 函数，在图形窗口显示十字光标，可以求得特殊点对应的K 值，进而分析系统稳定性情况。

已知一控制系统，H （s ）=1，其开环传递函数为：)

2)(1()(++=

s s s k s G ，绘制系统的根轨迹图，分析系统的稳定性。

程序为：

根轨迹图如图3所示：

图3 系统的根轨迹图光标选定虚轴临界点，程序结果为：

光标选定分离点，程序结果为：

上述数据显示了增益及对应的闭环极点位置。由此可得出如下结论：

（1）0

（2）k=0.4时，对应为分离点，系统处于临界阻尼状态；

（3）0.4

（4）k=6时，系统有一对虚根，系统处于临界稳定状态；

（5）k>6时，系统的一对复根的实部为正，系统处于不稳定状态。

3.4 Nyquist 曲线判断系统的稳定性

已知一控制系统，H （s ）=1，其开环传递函数为：)

2)(1()(++=s s s k s G 选用不同的增益K 值，用Matlab 绘出系统的阶跃响应曲线和Nyquist 曲线，并分析系统稳定性。

K=3，Nyquist 曲线程序：

K=9，Nyquist 曲线程序：

图4（a）K=3时Nyquist曲线

图4（a）K=9时Nyquist曲线

奈氏稳定判据的内容是：若开环传递函数在s平面右半平面上有P个极点，则当系统角频率X由-∞变到+∞时，如果开环频率特性的轨迹在复平面上逆时

针围绕（-1，j0）点转P圈，则闭环系统是稳定的，否则是不稳定的。

当k=3时，从图3（a）中可以看出，Nyquist曲线不包围（-1，j0）点，同时开环系统所有极点都位于s平面左半平面，因此，根据奈氏判据判定以此构成的闭环系统是稳定的；

当k=9时，从图3（b）中可以看出，Nyquist曲线按逆时针包围（-1，j0）点2圈，但此时P=0，所以根据奈氏判据判定以此构成的闭环系统是不稳定的。

4、结论

本文利用MATLAB函数的方法实现了对线性系统的稳定性分析,其过程简单,方法有效,结论直观,由此可见,MATLAB为工程技术人员分析、设计较优的控制系统提供了强有力的工具。现今，MATLAB也已成为大学里的许多先进课程的标准的计算机教学工具。相信不久将来,MATLAB也会出现在我国的许多大学和研究所里,并为各行各业专业人员所喜欢。

参考文献：

[1]、胡寿松.自动控制原理（第五版）[M].北京：科学出版社.2007.

[2]、张威.MATLAB基础与编程入门（第二版）[I].西安：西安电子科技大学出版社.2007.

[3]、高西全，丁美玉.数字信号处理（第三版）[I].西安：西安电子科技大学出版社.2008.

利用matlab进行系统分析基础

实验一利用matlab进行系统分析基础1．描述线性系统的三种不同方式之间的转换

问题1 已知系统的传递函数为 将其转换为零极点型。 相应的matlab语句为： num=[2 10]; den=[1 8 19 12]; printsys(num,den,’s’) 回车 [z,p,k]=tf2zp(num,den) 回车 察看语句的执行结果，并说明最后一行程序执行结果的含义；问题2 已知传递函数同上，试将其转换为状态变量型。Matlab语句为： Num=[2 10]; den=[1 8 19 12]; [a,b,c,d]=tf2ss(num,den) 回车

对应的状态方程为 式中A，B，C，D对应于程序中的a，b，c，d。

问题3 已知系统的零极点型传递函数为，试将其转换为传递函数型。Matlab语句： z=-1;p=[-2 –3 –4 ]; k=5; (回车) [num, den]=zp2tf(-1, [-2 –3 –4 ],2) （回车） %观察显示结果 继续输入： printsys(num,den,’s’) (回车) 记录显示结果。 2．卷积计算 原理： 两个信号卷积公式：

对于两个不规则波形的卷积，依靠手算是很困难的，在Matlab种则变得十分简单。 例如已知两个信号 其中分别表示两个门函数。 求其卷积的matlab程序如下： t1=1:0.01:2; f1=ones(size(t1)).*(t1>1);（表示一个高度为1的门函数，时间从t=1到 t=2） t2=2:0.01:3; f2=ones(size(t2)).*(t2>2); （表示一个高度为1的门函数，时间从t=2到t=3） c=conv(f1,f2);（卷积） t3=3:0.01:5; subplot(3,1,1),plot(t1,f1); subplot(3,1,2),plot(t2,f2); subplot(3,1,3),plot(t3,c); 其结果如图所示 问题1 已知两个信号 试利用matlab计算卷积 (要求显示出波形图） 3．傅立叶变换

用MATLAB实现线性系统的频域分析报告

实验二用MATLA实现线性系统的频域分析 [ 实验目的] 1 .掌握MATLAE平台下绘制典型环节及系统开环传递函数的Bode图和Nyquist图(极坐标图)绘制 方法； 2.掌握利用Bode图和Nyquist图对系统性能进行分析的理论和方法。 [ 实验指导] 一、绘制Bode图和Nyquist图 1.Bode图绘制 采用bode() 函数，调用格式： ①bode(sys) ; bode(num,den)； 系统自动地选择一个合适的频率围。 ②bode(sys , w)； 其中w(即3)是需要人工给出频率围，一般由语句w=logspace(a,b,n)给出。logspace(a,b,n):表示在10a到10b之间的n个点,得到对数等分的w值。 ③bode(sys,{wmin,wmax}); 其中{wmi n,wmax}是在命令中直接给定的频率w的区间。 以上这两种格式可直接画出规化的图形。 ④[mag,phase, 3 ]=bode(sys)或[m,p]=bode(sys) 这种格式只计算Bode图的幅值向量和相位向量，不画出图形。 m为频率特性G(j 3 )的幅值向量； p 为频率特性G(j 3 ) 的幅角向量，单位为角度(°)。 w为频率向量，单位为[弧度]/秒。在此基础上再画图，可用: subplot(211)；semilogx(w,20*log10(m) % 对数幅频曲线subplot(212)；semilogx(w,p) % 对数相频曲线 ⑤bode(sys1,sys2 ,…,sys N); ⑥bode((sys1,sys2 ,…,sys N, w); 这两种格式可在一个图形窗口同时绘多个系统的bode图。 2.Nyquist 曲线的绘制

电力电子技术与电力系统分析matlab仿真

电气2013级卓班电力电子技术与电力系统分析 课程实训报告 专业：电气工程及其自动化 班级： 姓名： 学号： 指导教师：

兰州交通大学自动化与电气工程学院 2016 年 1 月日

电力电子技术与电力系统分析课程实训报告 1 电力电子技术实训报告 1.1 实训题目 1.1.1电力电子技术实训题目一 一．单相半波整流 参考电力电子技术指导书中实验三负载，建立MATLAB/Simulink环境下三相半波整流电路和三相半波有源逆变电路的仿真模型。仿真参数设置如下： (1)交流电压源的参数设置和以前实验相关的参数一样。 (2)晶闸管的参数设置如下： R=0.001Ω，L =0H，V f=0.8V，R s=500Ω，C s=250e-9F on (3)负载的参数设置 RLC串联环节中的R对应R d，L对应L d，其负载根据类型不同做不同的调整。 (4)完成以下任务： ①仿真绘出电阻性负载（RLC串联负载环节中的R d= Ω，电感L d=0，C=inf，反电动势为0）下α=30°，60°，90°，120°，150°时整流电压U d，负载电流L 和晶闸管两端电压U vt1的波形。 d ②仿真绘出阻感性负载下（负载R d=Ω，电感L d为，反电动势E=0）α=30°，60°，90°，120°，150°时整流电压U d，负载电流L d和晶闸管两端电压U vt1的波形。 ③仿真绘出阻感性反电动势负载下α=90°，120°，150°时整流电压U d，负载电流L d和晶闸管两端电压U vt1的波形，注意反电动势E的极性。 (5)结合仿真结果回答以下问题： ①该三项半波可控整流电路在β=60°，90°时输出的电压有何差异?

实验一--控制系统的稳定性分析

实验一--控制系统的稳定性分析

实验一控制系统的稳定性分 班级：光伏2班 姓名：王永强 学号：1200309067

实验一控制系统的稳定性分析 一、实验目的 1、研究高阶系统的稳定性，验证稳定判据的正确性； 2、了解系统增益变化对系统稳定性的影响；

3、观察系统结构和稳态误差之间的关系。 二、实验任务 1、稳定性分析 欲判断系统的稳定性，只要求出系统的闭环极点即可，而系统的闭环极点就是闭环传递函数的分母多项式的根，可以利用MATLAB中的tf2zp函数求出系统的零极点，或者利用root函数求分母多项式的根来确定系统的闭环极点，从而判断系统的稳定性。 （1）已知单位负反馈控制系统的开环传递 函数为 0.2( 2.5) () (0.5)(0.7)(3) s G s s s s s + = +++，用MATLAB编写 程序来判断闭环系统的稳定性，并绘制闭环系统的零极点图。 在MATLAB命令窗口写入程序代码如下：z=-2.5 p=[0,-0.5,-0.7,-3] k=1 Go=zpk(z,p,k)

Gc=feedback(Go,1) Gctf=tf(Gc) dc=Gctf.den dens=ploy2str(dc{1},'s') 运行结果如下： Gctf = s + 2.5 --------------------------------------- s^4 + 4.2 s^3 + 3.95 s^2 + 2.05 s + 2.5 Continuous-time transfer function. dens是系统的特征多项式，接着输入如下MATLAB程序代码： den=[1,4.2,3.95,1.25,0.5] p=roots(den)

利用MATLAB进行时域分析

自动控制原理与系统课程实验报告 实验题目：利用MATLAB进行时域分析 班级：机电1131班姓名：刘润学号：38号 一、实验目的及内容 时域分析法是一种直接在时间域中对系统进行分析的方法，具有直观、准确的优点，并且可以提供系统时间响应的全部信息。在此实验中，主要介绍时域法进行系统分析，包括一阶系统、二阶系统以及高阶系统,以及系统的性能指标。通过实验，能够快速掌握、并利用MATLAB及控制系统箱对各种复杂控制系统进行时域分析。 二、实验设备 三、实验原理 典型的二阶系统在不同的阻尼比的情况下，它们的阶跃响应输出特性的差异是很大的。若阻尼比过小，则系统的振荡加剧，超调量大幅度增加；若阻尼比过大，则系统的响应过慢，又大大增加了调整时间，下面通过此实验课题分析输出响应变化规律： 已知二阶振荡环节的传递函数为：G（s）=ωn*ωn/（s*s+2*ζ*ωn*s+ωn*ωn）， 其中ωn=0.4，ζ从0变化到2，求此系统的单位阶跃响应曲线，并分析当ζ发生变化时，二阶系统的响应有什么样的变化规律。

四、实验步骤编出程序如下图： 五、实验结果画出图表如下图：

六、结果分析 （1）当ξ=0（无阻尼）（零阻尼）时： 无阻尼时的阶跃响应为等幅振荡曲线。如图ξ=0曲线。 （2）当0<ξ<1（欠阻尼）时： 对应不同的ξ，可画出一系列阻尼振荡曲线，且ξ越小，振荡的最大振幅愈大。如图ξ=0.4曲线。 (3)当ξ=1（临界阻尼）时： 临界阻尼时的阶跃响应为单调上升曲线。如图ξ=1曲线。 （4）当ξ>1（过阻尼）时： 过阻尼时的阶跃响应也为单调上升曲线。不过其上升的斜率较临界阻尼更慢。如图ξ=1.6曲线 七、教师评语

用Matlab计算潮流计算电力系统分析

《电力系统潮流上机》课程设计报告 院系：电气工程学院 班级：电088班 学号： 0812002221 学生姓名：刘东昇 指导教师：张新松 设计周数：两周 日期：2010年 12 月 25 日

一、课程设计的目的与要求 目的：培养学生的电力系统潮流计算机编程能力，掌握计算机潮流计算的相关知识 要求：基本要求： 1.编写潮流计算程序； 2.在计算机上调试通过； 3.运行程序并计算出正确结果； 4.写出课程设计报告 二、设计步骤： 1.根据给定的参数或工程具体要求（如图），收集和查阅资料；学习相关软件（软件自选：本设计选择Matlab进行设计）。 2.在给定的电力网络上画出等值电路图。 3.运用计算机进行潮流计算。 4.编写设计说明书。 三、设计原理 1．牛顿-拉夫逊原理 牛顿迭代法是取x0 之后，在这个基础上，找到比x0 更接近的方程的跟，一步一步迭代，从而找到更接近方程根的近似跟。牛顿迭代法是求方程根的重要方法之一，其最大优点是在方程f(x) = 0 的单根附近具有平方收敛，而且该法还可以用来求方程的重根、复根。电力系统潮流计算，一般来说，各个母线所供负荷的功率是已知的，各个节点电压是未知的（平衡节点外）可以根据网络结构形成节点导纳矩阵，然后由节点导纳矩阵列写功率方程，由于功率方程里功率是已知的，电压的幅值和相角是未知的，这样潮流计算的问题就转化为求解非线性方程组的问题了。为了便于用迭代法解方程组，需要将上述功率方程改写成功率平衡方程，并对功率平衡方程求偏导，得出对应的雅可比矩阵，给未知节点赋电压初值，一般为

额定电压，将初值带入功率平衡方程，得到功率不平衡量，这样由功率不平衡量、雅可比矩阵、节点电压不平衡量（未知的）构成了误差方程，解误差方程，得到节点电压不平衡量，节点电压加上节点电压不平衡量构成新的节点电压初值，将新的初值带入原来的功率平衡方程，并重新形成雅可比矩阵，然后计算新的电压不平衡量，这样不断迭代，不断修正，一般迭代三到五次就能收敛。 牛顿—拉夫逊迭代法的一般步骤： （1）形成各节点导纳矩阵Y。 （2）设个节点电压的初始值U和相角初始值e 还有迭代次数初值为0。 （3）计算各个节点的功率不平衡量。 （4）根据收敛条件判断是否满足，若不满足则向下进行。 （5）计算雅可比矩阵中的各元素。 （6）修正方程式个节点电压 （7）利用新值自第（3）步开始进入下一次迭代，直至达到精度退出循环。 （8）计算平衡节点输出功率和各线路功率 2．网络节点的优化 １）静态地按最少出线支路数编号 这种方法由称为静态优化法。在编号以前。首先统计电力网络个节点的出线支路数，然后，按出线支路数有少到多的节点顺序编号。当由n 个节点的出线支路相同时，则可以按任意次序对这n 个节点进行编号。这种编号方法的根据是导纳矩阵中，出线支路数最少的节点所对应的行中非零元素也２）动态地按增加出线支路数最少编号在上述的方法中，各节点的出线支路数是按原始网络统计出来的，在编号过程中认为固定不变的，事实上，在节点消去过程中，每消去一个节点以后，与该节点相连的各节点的出线支路数将发生变化(增加，减少或保持不变)。因此，如果每消去一个节点后，立即修正尚未编号节点的出线支路数，然后选其中支路数最少的一个节点进行编号，就可以预期得到更好的效果，动态按最少出线支路数编号方法的特点就是按出线最少原则编号时考虑了消去过程中各节点出线支路数目的变动情况。 3．MATLAB编程应用 Matlab 是“Matrix Laboratory”的缩写，主要包括：一般数值分析，矩阵运算、数字信号处理、建模、系统控制、优化和图形显示等应用程序。由于使用Matlab 编程运算与人进行科学计算的思路和表达方式完全一致，所以不像学习高级语言那样难于掌握，而且编程效率和计算效率极高，还可在计算机上直接输出结果和精美的图形拷贝，所以它的确为一高效的科研助手。 四、设计内容

matlab实验四 系统的零极点分析

实验四连续时间系统复频域分析和离散时间系统z域分析 一.实验目的： 1.掌握连续信号拉氏变换和拉氏反变换的基本实现方法。 2.熟悉laplace函数求拉普拉斯变换，ilaplace函数求拉氏反变换 的使用。 3.掌握用ztrans函数，iztrans函数求离散时间信号z变换和逆z 变换的基本实现方法。 4.掌握用freqs函数，freqz函数由连续时间系统和离散时间系统 系统函数求频率响应。 5.掌握zplane零极点绘图函数的使用并了解使用零极点图判断系 统稳定性的原理。 二、实验原理： 1.拉氏变换和逆变换 原函数()() ?象函数 f t F s 记作：[()]() =→拉氏变换 L f t F s 1[()]() -=→拉氏反变换 L F s f t 涉及函数：laplace,ilapace. 例如：

syms t;laplace(cos(2*t)) 结果为：ans =s/(s^2+4) syms s;ilaplace(1./(s+1)) 结果为：ans = exp(-t) 2. 系统传递函数H(s)或H(z)。 12121212...()()()...m m m n n n b s b s b B s H s A s a s a s a ----+++==+++ 112112...()()()...m m m n n n b z b z b B z H z A z a z a z a --+--++++==+++ 其中，B 为分子多项式系数，A 为分母多项式系数。 涉及函数：freqz,freqs. 3. 系统零极点分布与稳定性的判定。 对于连续时间系统，系统极点位于s 域左半平面，系统稳定。 对于离散时间系统，系统极点位于z 域单位圆内部，系统稳定。 涉及函数：zplane. 三、 实验内容 1． 验证性实验 a) 系统零极点的求解和作图

基于Matlab计算程序的电力系统运行分析课程设计

课程设计 课程名称：电力系统分析 设计题目：基于Matlab计算程序地电力系统运行分析学院：电力工程学院 专业：电气工程自动化 年级： 学生姓名： 指导教师： 日期： 教务处制

目录 前言 (1) 第一章参数计算 (2) 一、目标电网接线图 (2) 二、电网模型地建立 (3) 第二章潮流计算 (6) 一．系统参数地设置 (6) 二．程序地调试 (7) 三、对运行结果地分析 (13) 第三章短路故障地分析计算 (15) 一、三相短路 (15) 二、不对称短路 (16) 三、由上面表对运行结果地分析及在短路中地一些问题 (21) 心得体会 (26) 参考文献 (27)

前言 电力系统潮流计算是电力系统分析中地一种最基本地计算，是对复杂电力系统正常和故障条件下稳态运行状态地计算.潮流计算地目标是求取电力系统在给定运行状态地计算.即节点电压和功率分布，用以检查系统各元件是否过负荷.各点电压是否满足要求，功率地分布和分配是否合理以及功率损耗等.对现有电力系统地运行和扩建，对新地电力系统进行规划设计以及对电力系统进行静态和暂态稳定分析都是以潮流计算为基础.潮流计算结果可用如电力系统稳态研究，安全估计或最优潮流等对潮流计算地模型和方法有直接影响. 在电力系统中可能发生地各种故障中，危害最大且发生概率较高地首推短路故障.产生短路故障地主要原因是电力设备绝缘损坏.短路故障分为三相短路、两相短路、单相接地短路及两相接地短路.其中三相短路时三相电流仍然对称，其余三类短路统成为不对称短路.短路故障大多数发生在架空输电线路.电力系统设计与运行时，要采取适当地措施降低短路故障地发生概率.短路计算可以为设备地选择提供原始数据.

利用matlab分析系统动态性能

利用matlab分析系统动态性能

控制系统的时域分析 一.系统阶跃响应的性能指标 表 1 系统性能指标 利用 matlab 程序求出各系统阶跃响应的性能指标及图像，如求原系统 1 的方程： num=1.05; den=conv([0.125,1],conv([0.5,1],[1,1,1])); G=tf(num,den); C=dcgain(G); [y,t]=step(G); plot(t,y) grid [Y,K]=max(y); tp=t(K) mp=100*(Y-C)/C n=1; while y(n)0.98*C)&&(y(i)<1.02*C) i=i-1; end ts=t(i)

图 1 系统 1 阶跃响应曲线图二．根据系统性能指标及图像分析系统 1.利用 Matlab 得各系统节约系统曲线，如图 2：num1=1.05; den1=conv([0.125,1],conv([0.5,1],[1,1,1])); G1=tf(num1,den1); [y1,t1]=step(G1); num2=1.05*[0.4762,1]; den2=conv([0.125,1],conv([0.5,1],[1,1,1])); G2=tf(num2,den2); [y2,t2]=step(G2); num3=1.05*[1,1]; den3=conv([0.125,1],conv([0.5,1],[1,1,1])); G3=tf(num3,den3); [y3,t3]=step(G3); num4=1.05*[0.4762,1]; den4=conv([0.25,1],conv([0.5,1],[1,1,1])); G4=tf(num4,den4); [y4,t4]=step(G4); num5=1.05*[0.4762,1]; den5=conv([0.5,1],[1,1,1]); G5=tf(num5,den5); [y5,t5]=step(G5); num6=1.05; den6=[1,1,1]; G6=tf(num6,den6);

分析非线性系统的方法

非线性系统稳定性问题的判定方法和发展趋势 任何一个实际系统总是在各种偶然和持续的干扰下运动或工作的。所以，当系统承受干扰之后，能否稳妥地保持预订的运动轨迹或者工作状态，即系统的稳定性是首要考虑的。一个系统的稳定性，包括平衡态的稳定性问题和任一运动的稳定性问题。而对于给定运动的稳定性可以变换成关于平衡点的稳定性问题。 对平衡点的稳定性进行分析可将平衡点的稳定性定义为李雅普诺夫稳定、一致稳定、渐进稳定、一致渐近稳定、按指数渐进稳定和全局渐进稳定，除了全局渐进稳定，其他都是局部的概念。 非线性系统的数学模型不满足叠加原理或其中包含非线性环节。包括非本质非线性（能够用小偏差线性化方法进行线性化处理的非线性）和本质非线性（用小偏差线性化方法不能解决的非线性）。它与线性系统有以下主要区别： 1．线性控制系统只能有一个平衡点或无穷多的平衡点。但非线性系统可以有一个、二个、多个、以至无穷多个平衡点。非线性系统与线性定常系统明显不同，其稳定性是针对各个平衡点而言的。通常不能说系统的稳定性如何，而应说那个平衡点是稳定的或不稳定的。2．在线性系统中，系统的稳定性只与系统的结构和参数有关，而与外作用及初始条件无关。非线性系统的稳定性除了与系统的结构和参数有关外，还与外作用及初始条件有关。 由于非线性控制系统与线性控制系统有很大的差异，因此，不能直接用线性理论去分析它，否则会导致错误的结论。对非线性控制系统的分析，还没有一种象线性控制系统那么普遍的分析、设计方法。 现代广泛应用于非线性系统上的分析方法有基于频率域分析的描述函数法和波波夫超稳定性，还有基于时间域分析的相平面法和李雅普诺夫稳定性理论等。这些方法分别在一定的假设条件下，能提供关于系统稳定性或过渡过程的信息。而计算机技术的迅速发展为分析和设计复杂的非线性系统提供了有利的条件。另外，在工程上还经常遇到一类弱非线性系统，即特性和运动模式与线性系统相差很小的系统。对于这类系统通常以线性系统模型作为一阶近似，得出结果后再根据系统的弱非线性加以修正，以便得到较精确的结果。摄动方法是处理这类系统的常用工具。而对于本质非线性系统，则需要用分段线性化法等非线性理论和方法来处理。目前分析非线性控制系统的常用方法如下： 1、线性化方法 采用线性化模型来近似分析非线性系统。 这种近似一般只限于在工作点附近的小信号情况下才是正确的。这种线性化近似，只是对具有弱非线性（或称非本质非线性）的系统。 常用线性化方法，有正切近似法和最小二乘法。 此外，对一些物理系统的非线性特性比较显著，甚至在工作点附件的小范围内也是非线性的，并且不能用一条简单的直线来代表整个非线性系统特性的系统，可采用分段线性化方法。2、相平面法 相平面法是一种基于时域的分析方法，一种用图解法求解一、二阶非线性常微分方程的方法。 该方法通过图解法将一阶和二阶系统的运动过程转化为位置和速度平面上的相轨迹，从而比较直观、准确地反映系统的稳定性、平衡状态和稳态精度以及初始条件及参数对系统运动的影响。相轨迹的绘制方法步骤简单、计算量小，特别适用于分析常见非线性特性和一阶、二阶线性环节组合而成的非线性系统 对于分段线性的非线性系统来说，相平面分析法的步骤为： （1）用n条分界线（开关线，转换线）将相平面分成n个线性区域；（2）分别写出各个线性区域的微分方程；（3）求出各线性区的奇点位置并画出相平面图；

非线性系统稳定性问题的判定方法和发展趋势

非线性系统的概念及稳定性问题的判定方法和发展趋势 姓名：查晓锐 学号：121306060006 线性系统理论自20世纪50年代以来不仅已在理论上逐步完善，也已成功的应用于各种国防和工业控制问题。随着现代工业对控制系统性能的要求不断提高，传统的线性反馈控制已很难满足各种实际需要。这是因为大多数实际控制系统往往是非线性的，采用近似的线性模型虽然可以使我们更全面和容易的分析系统的各种特性，但是却很难刻画出系统的非线性本质，线性系统的动态特性已不足以解释许多常见的实际非线性现象。另一方面，计算机及传感器技术的飞速发展，也为我们实现各种复杂非线性控制算法奠定了硬件基础。因此自20世纪80年代以来，非线性系统的控制问题受到了国内外控制界的普遍关注。 非线性科学是当今世界科学的前沿与热点，涉及自然科学和人文社会科学的众多领域，具有重大的科学价值和深刻的哲学方法论意义。但迄今为止，对非线性的概念、非线性的性质，并没有清晰的、完整的认识，对其哲学意义也没有充分地开掘。 一、 非线性的概念 非线性是相对于线性而言的，对线性的否定，线性是非线性的特例。所以要弄清非线性的概念，明确什么是非线性，首先必须明确什么是线性；其次对非线性的界定必须从数学表述和物理意义两个方面阐述，才能较完整地理解非线性的概念。 对线性的界定，一般是从相互关联的两个角度来进行的。其一：叠加原理成立“ 如果1Φ，2Φ 是两个那么21Φ+Φβα也是它的一个解，换言之，两个态的叠加仍然是一个态。”原理成立意味着所考查系统的子系统间没有非线性相互作用。其二，物理变量间的函数关系是直线，变量间的变化率是恒量，这意味着函数的斜率在其定义域内处处存在且相等，量间的比例关系在变量的整个定义域内是对称的。 在明确了线性的含义后，相应地非线性概念就易于界定。其一 ：“定义非线性算符()ΦN 为对一些 a ，b 或Φ，ψ不满足)()()(ψ+Φ=ψ+ΦbL aL b a L 的算符 即叠加原理不成立。”这意味着Φ与ψ之间存在藕合，对ψ+Φb a 的操作，等于分别对Φ，ψ操作外，再加上对Φ与ψ的交叉项（耦合项）操作，或者Φ、ψ是不连续有突变或断裂、不可微有折点的。其二：作为等价的另一种表述，我们可以从另一个角度来理解非线性在用于描述一个系统的一套确定的物理变量中，一个系统的一个变量最初的变化所造成的此变量或其它变量的相应变化是不成比例的。换言之：变量间的变化率不是恒量，函数的斜率在其定义域中有不存在或不相等的地方。概括地说：物理变量间的一级增量关系在变量的定义域内是不

(完整版)电力系统分析大作业matlab三机九节点潮流计算报告

电力系统分析大作业

一、设计题目 本次设计题目选自课本第五章例5-8,美国西部联合电网WSCC系统的简化三机九节点系统，例题中已经给出了潮流结果，计算结果可以与之对照。取ε=0.00001 。

二、计算步骤 第一步，为了方便编程，修改节点的序号，将平衡节点放在最后。如下图： 第二步，这样得出的系统参数如下表所示： 第三步，形成节点导纳矩阵。 9 2 1 3 2 7 4 5 6 8 3

第四步，设定初值： ο 01)0(6)0(5)0(4)0(3)0(2)0(1∠======??????U U U U U U ； 0)0(8)0(7==Q Q ，0)0(8)0(7==θθ。 第五步，计算失配功率 )0(1P ?=0，)0(2P ?=-1.25，)0(3P ?=-0.9，) 0(4P ?=0，)0(5P ?=-1，)0(6P ?=0，)0(7P ?=1.63， )0(8P ?=0.85； )0(1Q ?=0.8614，)0(2Q ?=-0.2590，)0(3Q ?=-0.0420，) 0(4Q ?=0.6275，)0(5Q ?=-0.1710， )0(6Q ?=0.7101。 显然，5108614.0|},max {|-=>=??εi i Q P 。 第六步，形成雅克比矩阵（阶数为14×14） 第七步，解修正方程，得到： =?)0(1θ-0.0371，=?)0(2θ-0.0668，=?)0(3θ-0.0628，=?)0(4θ0.0732，=?)0(5θ0.0191，=?)0(6θ0.0422，=?)0(7θ0.1726，=?)0(8θ0.0908； =?)0(1U 0.0334，=?)0(2U 0.0084，=?)0(3U 0.0223，=?)0(4U 0.0372，=?)0(5U 0.0266，

用Matlab计算潮流计算电力系统分析

《电力系统潮流上机》课程设计报告院系：电气工程学院 班级：电088班 学号： 学生姓名：刘东昇 指导教师：张新松 设计周数：两周 日期：2010年 12 月 25 日

一、课程设计的目的与要求 目的：培养学生的电力系统潮流计算机编程能力，掌握计算机潮流计算的相关知识 要求：基本要求： 1.编写潮流计算程序； 2.在计算机上调试通过； 3.运行程序并计算出正确结果； 4.写出课程设计报告 二、设计步骤： 1.根据给定的参数或工程具体要求（如图），收集和查阅资料；学习相关软件（软件自选：本设计选择Matlab进行设计）。 2.在给定的电力网络上画出等值电路图。 3.运用计算机进行潮流计算。 4.编写设计说明书。 三、设计原理 1．牛顿-拉夫逊原理 牛顿迭代法是取x0 之后，在这个基础上，找到比x0 更接近的方程的跟，一步一步迭代，从而找到更接近方程根的近似跟。牛顿迭代法是求方程根的重要方法之一，其最大优点是在方程f(x) = 0 的单根附近具有平方收敛，而且该法还可以用来求方程的重根、复根。电力系统潮流计算，一般来说，各个母线所供负荷的功率是已知的，各个节点电压是未知的（平衡节点外）可以根据网络结构形成节点导纳矩阵，然后由节点导纳矩阵列写功率方程，由于功率方程里功率是已知的，电压的幅值和相角是未知的，这样潮流计算的问题就转化为求解非线性方程组的问题了。为了便于用迭代法

解方程组，需要将上述功率方程改写成功率平衡方程，并对功率平衡方程求偏导，得出对应的雅可比矩阵，给未知节点赋电压初值，一般为额定电压，将初值带入功率平衡方程，得到功率不平衡量，这样由功率不平衡量、雅可比矩阵、节点电压不平衡量（未知的）构成了误差方程，解误差方程，得到节点电压不平衡量，节点电压加上节点电压不平衡量构成新的节点电压初值，将新的初值带入原来的功率平衡方程，并重新形成雅可比矩阵，然后计算新的电压不平衡量，这样不断迭代，不断修正，一般迭代三到五次就能收敛。 牛顿—拉夫逊迭代法的一般步骤： （1）形成各节点导纳矩阵Y。 （2）设个节点电压的初始值U和相角初始值e 还有迭代次数初值为0。 （3）计算各个节点的功率不平衡量。 （4）根据收敛条件判断是否满足，若不满足则向下进行。 （5）计算雅可比矩阵中的各元素。 （6）修正方程式个节点电压 （7）利用新值自第（3）步开始进入下一次迭代，直至达到精度退出循环。 （8）计算平衡节点输出功率和各线路功率 2．网络节点的优化 １）静态地按最少出线支路数编号 这种方法由称为静态优化法。在编号以前。首先统计电力网络个节点的出线支路数，然后，按出线支路数有少到多的节点顺序编号。当由n 个节点的出线支路相同时，则可以按任意次序对这n 个节点进行编号。这种编号方法的根据是导纳矩阵中，出线支路数最少的节点所对应的行中非零元素也

控制系统的稳定性分析

精品 实验题目控制系统的稳定性分析 一、实验目的 1．观察系统的不稳定现象。 2．研究系统开环增益和时间常数对稳定性的影响。 二、实验仪器 1．EL-AT-II型自动控制系统实验箱一台 2．计算机一台 三、系统模拟电路图 系统模拟电路图如图3-1 图3-1 系统模拟电路图R3=0～500K； C=1μf或C=0.1μf两种情况。 四、实验报告 1．根据所示模拟电路图，求出系统的传递函数表达式。 G(S)= K=R3/100K,T=CuF/10 2．绘制EWB图和Simulink仿真图。

精品 3．根据表中数据绘制响应曲线。 4．计算系统的临界放大系数，确定此时R3的值，并记录响应曲线。 系统响应曲线 实验曲线Matlab （或EWB）仿真 R3=100K = C=1UF 临界 稳定 （理论值 R3= 200K） C=1UF

精品 临界 稳定 （实测值 R3= 220K） C=1UF R3 =100K C= 0.1UF

精品 临界 稳定 （理论 值R3= 1100 K） C=0.1UF 临界稳定 （实测值 R3= 1110K ） C= 0.1UF

精品 实验和仿真结果 1．根据表格中所给数据分别进行实验箱、EWB或Simulink实验，并进行实验曲线对比，分析实验箱的实验曲线与仿真曲线差异的原因。 对比： 实验曲线中R3取实验值时更接近等幅振荡，而MATLAB仿真时R3取理论值更接近等幅振荡。 原因： MATLAB仿真没有误差，而实验时存在误差。 2．通过实验箱测定系统临界稳定增益，并与理论值及其仿真结果进行比较（1）当C=1uf，R3=200K(理论值)时，临界稳态增益K=2， 当C=1uf，R3=220K(实验值)时，临界稳态增益K=2.2，与理论值相近（2）当C=0.1uf，R3=1100K(理论值)时，临界稳态增益K=11 当C=0.1uf，R3=1110K(实验值)时，临界稳态增益K=11.1，与理论值相近 四、实验总结与思考 1．实验中出现的问题及解决办法 问题：系统传递函数曲线出现截止失真。 解决方法：调节R3。 2．本次实验的不足与改进 遇到问题时，没有冷静分析。考虑问题不够全面，只想到是实验箱线路的问题，而只是分模块连接电路。 改进：在实验老师的指导下，我们发现是R3的取值出现了问题，并及时解决，后续问题能够做到举一反三。 3．本次实验的体会 遇到问题时应该冷静下来，全面地分析问题。遇到无法独立解决的问题，要及时请教老师，

利用MATLAB软件分析系统的频率响应

备注：（1）、按照要求独立完成实验项目内容，报告中要有程序代码和程序运行结果和波形图等原始截图。 （2）、实验结束后，把电子版实验报告按要求格式改名（例：09号-张三-实验一）后，上传至指定ftp服务器目录下（homework_upload）的相应文件里，并由实验教师批阅记录后； 实验室统一刻盘留档。 ftp:59.74.50.66 账号：microele 密码：ele1507 实验七利用MATLAB软件分析系统的频率响应 一、实验目的： 1、利用MATLAB求解系统的频率响应。 二、实验原理 MATLAB提供了函数freqs来计算系统的频率响应。 三、实验内容：（包括代码与产生的图形） 6-16 w=linspace(0,5,200); b=[1]; a=[1 2 2 1]; H=freqs(b,a,w); subplot(2,1,1); plot(w,abs(H)); set(gca,'xtick',[0 1 2 3 4 5]); set(gca,'ytick',[0 0.4 0.7071]);grid; xlabel('\omega') subplot(2,1,2); plot(w,angle(H)); set(gca,'xtick',[0 1 2 3 4 5]);grid; xlabel('\omega');

012345 00.4 0.7071 ω 012345 -4-2 24 ω 6-17 RC=0.04; t=linspace(-2,2,1024); w1=5;w2=100; H1=j*w1/(j*w1+1/RC); H2=j*w2/(j*w2+1/RC); f=cos(5*t)+cos(100*t); y=abs(H1)*cos(w1*t+angle(H1))+abs(H2)*cos(w2*t+angle(H2)); subplot(2,1,1); plot(t,f); subplot(2,1,2); plot(t,y);

实验五 线性系统的稳定性和稳态误差分析

实验五 自动控制系统的稳定性和稳态误差分析 一、实验目的 1、研究高阶系统的稳定性，验证稳定判据的正确性； 2、了解系统增益变化对系统稳定性的影响； 3、观察系统结构和稳态误差之间的关系。 二、实验任务 1、稳定性分析 欲判断系统的稳定性，只要求出系统的闭环极点即可，而系统的闭环极点就是闭环传递函数的分母多项式的根，可以利用MATLAB 中的tf2zp 函数求出系统的零极点，或者利用root 函数求分母多项式的根来确定系统的闭环极点，从而判断系统的稳定性。 （1）已知单位负反馈控制系统的开环传递函数为 0.2( 2.5) ()(0.5)(0.7)(3) s G s s s s s += +++，用MATLAB 编写程序来判断闭环系统的稳定性， 并绘制闭环系统的零极点图。 在MATLAB 命令窗口写入程序代码如下： z=-2.5 p=[0,-0.5,-0.7,-3] k=0.2 Go=zpk(z,p,k) Gc=feedback(Go,1) Gctf=tf(Gc) dc=Gctf.den dens=poly2str(dc{1},'s') 运行结果如下： dens= s^4 + 4.2 s^3 + 3.95 s^2 + 1.25 s + 0.5 dens 是系统的特征多项式，接着输入如下MATLAB 程序代码： den=[1,4.2,3.95,1.25,0.5]

p=roots(den) 运行结果如下： p = -3.0058 -1.0000 -0.0971 + 0.3961i -0.0971 - 0.3961i p为特征多项式dens的根，即为系统的闭环极点，所有闭环极点都是负的实部，因此闭环系统是稳定的。 下面绘制系统的零极点图，MATLAB程序代码如下： z=-2.5 p=[0,-0.5,-0.7,-3] k=0.2 Go=zpk(z,p,k) Gc=feedback(Go,1) Gctf=tf(Gc) [z,p,k]=zpkdata(Gctf,'v') pzmap(Gctf) grid 运行结果如下： z = -2.5000 p = -3.0058 -1.0000 -0.0971 + 0.3961i -0.0971 - 0.3961i k = 0.2000

利用MATLAB仿真软件系统进行图像的数据分析

课程设计任务书 学生姓名：叶枫专业班级：通信zy1201班指导教师：姜宁工作单位：信息工程学院 题目: 利用MATLAB仿真软件系统进行图像的数据分析初始条件： 1.MATLAB软件。 2.数字信号处理与图像处理基础知识。 要求完成的主要任务: 读取图像并求出图像的最大值、最小值、均值、中值、和、标准差、两图像的协方差、相关系数等。 课程设计的目的： 1．理论目的 课程设计的目的之一是为了巩固课堂理论学习，并能用所学理论知识正确分析信号处理的基本问题和解释信号处理的基本现象。 2．实践目的 课程设计的目的之二是通过设计具体的图像信号变换掌握图像和信号处理的方法和步骤。 时间安排： 指导教师签名：年月日 系主任（或责任教师）签名：年月日

目录 摘要 (3) Abstract (4) 1.Matlab及课程设计所用函数简介 (5) 1.1Matlab简介 (5) 1.2课程设计所用函数简介 (6) 2.数据采集 (9) 2.1 MATLAB的读取方法 (9) 3图像数据统计处理 (13) 3.1 图像数据处理原理 (13) 3.2各像素点中最大值的获取 (14) 3.3各像素点中最小值的获取 (14) 3.4各像素点值的均值的获取 (15) 3.5各像素点值的中值的获取 (16) 3.6各像素点值的和的获取 (17) 3.7各像素点值的标准差的获取 (18) 3.8各像素点值的方差的获取 (19) 3.9两图中各像素点值的协方差的获取 (20) 3.10两图的相对系数的获取 (20) 4.心得体会 (22) 参考文献 (23)

摘要 MATLAB软件是矩阵实验室的简称，是美国MathWorks公司出品的商业数学软件，可用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，广泛用于数字信号分析，系统识别，时序分析与建模，神经网络、动态仿真等方面有着广泛的应用。 MATLAB因具有强大的图形处理功能、符号运算功能和数值计算功能，而被广泛应用。而且随着信息时代和数字世界的到来，数字信号处理也已成为当今一门极其重要的学科和技术领域。目前数字信号处理在通信、语音、图像、自动控制、雷达、军事、航空航天、医疗和家用电器等众多领域得到了广泛的应用。 本次课程设计利用MATLAB软件进行图像的数据分析，包括读取图像并求出图像的最大值、最小值、均值、中值、和、标准差、两图像的协方差、相关系数等。 关键词：MATLAB，数字信号处理，图像数据分析

MATLAB经典控制系统的分析和设计

9.3.1 基本操作及命令 1. 访问和退出MATLAB 在大多数系统中，一旦安装了MATLAB，在调用时，应执行命令MATLAB。退出MATLAB应执行命令exit或quit。 2. 如何应用MATLAB 通常以命令驱动方式应用MATLAB。当输入单个命令时，MATLAB会立即对其进行处理，并且显示处理结果。MATLAB 也能够执行存储在文件中的命令序列。 通过键盘输入的命令，应用向上箭头键可以被存取。通过输入某个最新命令和调用特定的命令行，可以使屏幕内容向上滚动。 3. MATLAB的变量 MATLAB的一个特点是变量在应用之前不必是维数确定的。在MATLAB中，变量一旦被采用，会自动产生（如果必要，变量的维数以后还可以改变）。在命令exit或quit输入之前，这些变量将保留在存储器中。 为了得到工作空间内的变量清单，可以通过键盘输入命令who，当前存放在工作空间内的所有变量便会显示在屏幕上。 命令clear能从工作空间中清除所有非永久性变量。如果只需要从工作空间中清除某个特定变量，比如“x”，则应输入命令clear x。 4. 以“%”开始的程序行 在本书中，许多MATLAB程序在编写时附有注解和说明，这些注解和说明阐明了发生在程序中的具体进程。在MATLAB 中以“%”开始的程序行，表示注解和说明。符号“%”类似于BASIC中“REM”。以“%”开始的行，用来存储程序的注解或说明，这些注解和说明是不执行的。这就是说，在MATLAB程序行中，出现在“%”以后的一切内容都是可以忽略的。如果注解或说明需要一行以上程序行，则每一行均需以“%”为起始。 5. 应用分号操作符 分号用来取消打印。如果语句的最后一个符号是分号，则打印被取消，但是命令仍在执行，而结果不再显示。这是一个有益的特性，因为打印中间结果可能不必要。此外，在输入矩阵时，除非最后一行，分号用来指示一行的结束。 6. 应用冒号操作符

电力系统稳定性分析matlab程序

电力系统稳定性分析作业一 1 euler.m ,reuler.m, kunta.m分别为（1）中的欧拉法，改进欧拉法，龙格库塔法的主程序；doty.m,doty2.m,doty3.m均为（1）中子函数程序。Runge-Kutta.m为（2）和（3）的运行程序。 下表为三种方法的部分运行结果功角数据： 时 间 00.010.020.030.040.050.060.070.08 改进35.161 5 35.222 2 35.402 3 35.699 9 36.113 36.639 4 37.277 38.023 4 38.876 3 时 间 0．090.100.110.120.130.140.150.160.17 改39.83340．8942.00543.12544.25045.37546.49947.61848.730

进 2 18 1 6 1 7 5 7 5 （1）欧拉法 在matlaB中输入命令[t,x,y,z]=euler('doty','doty2','doty3',0,5,0.1,0.01)可得 t-w曲线，t-δ曲线分别如下图所示。具体功角，角速度数据分别见文件1.mat 和2.mat （2）欧拉改进法 在matlab命令窗口输入[t,x,y,z]=reuler('doty','doty2','doty3',0,5,0.1,0.01) t-w曲线，t-δ曲线分别如下图所示。具体功角，角速度数据分别见文件3.mat 和4.mat

（3）龙格库塔法 在matlab命令窗口输入[t,x,y,z]=kunta('doty','doty2','doty3',0,5,0.1,0.01) t-w曲线，t-δ曲线分别如下图所示。具体功角，角速度数据分别见文件5.mat 和6.mat 2 运行Runge-Kutta,将参数阻尼D设置为0.05，不断更改参数切除时间t的值，当t=0.2728和t=0.2730时，运行程序分别得到如下两图：