第8章电机的状态变量分析法

1、第八章 电机的状态变量分析法本世纪二十年代以来，研究连续信号作用于线性网络问题的最有力工具是拉氏变换，其特点是将研究对象的特性用传递函数表示，然后根据激励函数求出响应函数。这种方法称为输入词输出法，又称端部法，因为它只研究网络的端口特性，而不考虑网络内部的结构、参数、电压、电流等。六十年代以来，随着电子计算机的广泛应用，状态变量分析法有了很大的发展。与输入输出法相比，这种方法是一种内部法，因为它首先选定能够代表网络(或物理系统)内部特性的某些物理量，作为状态变量，建立相应的状态方程求解，然后由解出的状态变量和给定的输入量求出所需的输出量。机电系统(包括电机)可化成电网络求解。使用状态变量法来分

2、析机电系统的动态性能有以下优点：(1)状态变量方程是一组联立的一阶微分方程，可不必经过拉氏变换而直接在时域内求解，用电子计算机解这类方程组早已规格化，用起来很方便。(2)不仅适用于单输入、单输出的问题，也适用于多输入、多输出的问题，后者用输入呻输出法较难求解。(3)对于较复杂的网络，包括一部分线性、参数不随时间变化的网络和大部分非线性、时变参数的网络，数值计算法常是唯一可采用的方法。在用电子计算机进行这种计算时，用状态变量法要比用拉氏变换法容易。(4)由于状态变量法是一种内部法，用于分析系统的稳定性、可控度等是很方便的，可以在不求出变量的解的情况下进行。由于上述优点，用状态变量法分析电机的动态

3、性能，不仅有可能建立较精确的数学模型，使计算结果接近于实测值；而且可以解决某些过去难以分析的动态课题，从而使电机理论得到进一步发展。本章首先简要介绍状态变量和状态方程的基本概念，然后建立几种典型的电机的状态方程，最后概略介绍它们的求解方法，并附有计算实例。81状态变量分析法的基本概念一、状态变量和状态方程的基本概念由电路理论可知电容和电感是储能元件，它在某一瞬时的输出量不仅决定于该瞬时的输入，而是要由输入的全部历史来决定。例如电容器的电荷(输出量)与电流(输入量)的关系是 （81）这里是输入开始作用的时刻。在此之前，电容器中无电荷储存，即当，。由上式可见，为了确定在某一瞬时的输出值，需要知道在

4、区间的全部输入值。因此储能元件又可称为记忆元件，而非储能元件可称为无记忆元件。但是事实上不必从输入开始作用的时间分析起。对于任一给定时刻 (这时输入已开始了一段时间)，只要知道了该瞬时电荷的数值，即初始值，即可计算时的电荷 （82）而并不需要知道之前的和变化情况。就是含有电容的一阶网络的状态变量，是该状态变量的初始值，由它和在以后的外施激励，即可确定后网络的变化过程。网络中的状态变量可以不止一个。一般说来对应每一储能元件可选定一状态变量。在建立状态方程时，应使状态变量的数目为最小。因此，网络的状态变量就是这样一组描述该网络状态所须具备的最少变量，只要知道在某一初始时刻时的状态变量值和在区间网络

5、的输入量就可以完全确定在区间该网络的状态变量值和所要求出的输出量 ， 换句话说，状态变量概括了为预测网络未来的特性而必须知道的有关网络过去状况的信息。下面用一简单实例复习电路课程中学过的状态方程的概念。图8l为一RLC电路，R、L、C为已知参数，和为输入量，为输出量。 图81 RLC电路此电路的回路方程为 （83）这是一个二阶微分方程。用来确定积分常数的初始条件，是开关合下时刻时电容上电压和电感中电流。如改用电容电压和电感电流为变量来列上述电路的方程，则有 将这两个方程改写，使每个状态变量的一阶时间导数处于等式左方，而右方只含状态变量和输入变量而不含其导数，可得 （84）写成矩阵形式为 （85

6、）上式是一组以和为变量的联立的一阶微分方程，变量的初始值为和；所以和就是上述电路的状态变量，而一阶微分方程组（85）就是描述电路变化过程的状态方程。输出变量与状态变量之间的关系式，称为输出方程，在本例中 （86）即为输出方程。普遍地说，线性系统的状态方程和输出方程可写成如下标准形式 （87）式中为状态向量，它是由个状态变量所组成的列向量；为输入向量（或称控制向量），它是由个输入变量所组成的列向量；为输出向量，它是由个输出变量所组成的列向量；为的系统矩阵；为的控制矩阵；为的输出矩阵；为的矩阵。对于此电路来说，,从式（8-5）和（8-6）可见 （88） （89）如果网络中含有非线性元件，则其特性不

7、能简单地用常数、来代表，而只能表示为某种非线性函数，如等这时各状态变量的导数，不再等于状态变量和输入变量的线性组合，而只能表示为它们的某种非线性函数。因此对于非线性系统，状态方程和输出方程的一般形式为 （810）式中是一组个非线性函数，是一组个非线性函数。最后还应指出，状态变量的选择不是唯一的，例如储能元件可以取其磁链或其电流或其电压作为状态变量。一般来说，对于线性电路系统，可选电容电压和电感电流作为状态变量；对于非线性电路系统，则常选电容电荷和电感磁链作为状态变量。对于机械系统，一般选择角位移（位移）和转速（速度）作为状态变量。二、状态方程的建立步骤对于简单的系统，可按下述步骤建立状态方程：

8、（1） 根据系统情况选择合适的状态变量；（2） 运用基尔霍夫定律、大朗贝尔原理（或用拉格朗日方程）列写系统的运动方程，并进行适当简化；（3） 当方程中含有非状态变量时，应设法将其消去；（4） 将状态变量的一阶导数项移到方程的左方，将状态变量和输入变量项移到右方，将方程整理成状态方程的标准形式。有时已给出所描述研究对象的高阶微分方程，可设法将其转化成对应的一组状态方程。例如已知三阶微分方程为 （811）可改写为在上式中，令，即可化出一组状态方程为 （812a）也就是将对的一阶导数定义为第二个状态变量，将其对的二阶导数定义为第三个状态变量。上面的方程组写成矩阵形式时，则为 （812b）这就是三阶微

9、分方程（8-11）改写出的状态方程的标准形式。对于更为复杂的系统，可按网络图论来建立状态方程（有关内容可参阅专门文献，对于电机一般不需要）。8-2 电机的状态方程下面分别对直流电机、感应电机和同步电机的状态方程作一扼要介绍。一、 直流电机的状态方程以图8-2所示他励直流电动机为例说明如何建立状态方程。该电机的输入量为电枢电压、励磁电压和负载转矩，输出量为角位移和角速度。图82 他励直流电动机选择电枢电流、励磁电流和角速度作为状态变量。直流电动机的运动方程为 （8-13）式中为电枢感应电动势，为电磁转矩，为电机转子和负载的转动惯量，为旋转阻力系数。为了消去方程组中的非状态变量和，根据直流电机基本

10、理论有 （814）将式（812）代入式（815），得 （815）将上式整理后，可得 （816）写成矩阵形式为 （817）其中 （818）上式是直流电动机标准形式的状态方程。由于方程式右方有两个状态变量的乘积项，使得系数矩阵A中也含有状态变量，所以它是一组非线性微分方程，通常只能采用数值积分法求解。二、 感应电动机的状态方程以三项绕线式感应电动机为例。一般选择定、转子电流，角速度和电角位移等八个变量作为状态变量。各物理量的正方向采用电动机惯例，列写运动方程为 （819）式中为电机极对数，为电磁转矩。展开上式可得 （820） （821）式中下标和分别表示定子和转子，表示矩阵的转置。当定子三相绕组接

11、上三相对称正弦电压，而转子三相绕组自行短接时，式中各子矩阵分别为 （822） （823） (824) 由于感应电机气隙均匀，在忽略定、转子表面齿槽影响和铁心饱和的情况下，可假定气隙磁场按正弦分布，则定子和转子绕组的自感以及定子或转子绕组间的互感均为常量，与转子位置无关。绕组的轴线在空间互差120°电角度， 为负，使互感为负值，故得 (825) (826) 这样式（820）和（821）中定子和转子的电感子矩阵分别为 （827）由于定转子之间有相对运动，因而定子任一绕组和转子任一绕组间的互感将随该两绕组轴线间夹角而改变。当定、转子两个相绕组轴线重合时，他们之间的互感达到最大值。假定转子以

12、转速逆时钟向旋转且轴和轴之间的电弧度为，则 （828）式中为时的初始值。由于假定感应电机的气隙磁通密度在空间按正弦分布，这时定、转子两个相绕组间的互感将随该两绕组轴线间夹角的余弦函数而变化，于是式 （820）和（821）中的定、转子的互感矩阵为 （829）根据式（821），并考虑到，则电磁转矩应为 （830）式（819）是一组变系数的非线性微分方程，需要用数值法求解，可将其改写为 （831）即为标准形式的状态方程 （8-32）式中 （833） （834）上面是感应电机用abc坐标系统列写的状态方程，为简化计算，也可采用综合矢量表示。由式（567）可写出感应电机电压综合矢量方程 (835) 上式

13、可改写为 （836）或简作 即 （837）感应电机的转矩方程由式（558）为 若以电角速度为状态变量则上式可改写为 （838）式中由式（584）可写成 （839）将式（837）和（838）合并整理后可得状态变量、和表示的感应电机的运动方程为 （840）上式中多数项是线性的，也有一部分是非线性的，将其写成标准形式为式中而是非线性项。三、同步电机的状态同步电机可按和等坐标系统，采用磁链或电流作为状态变量来写状态方程。本节按较常采用的坐标系统、电动机惯例和基值系统并取电流为状态变量来研究这一问题。假定电机中不存在零序电流，则由式（631）、（632）和（633）可写出同步电机的电压和磁链方程如下 （

14、841）和 （842）利用上式及其两侧对时间求导数后的关系代入式（841）得 （843）上式也可改写为 （844）由式（636）同步电机的转矩方程为如取电角速度为状态变量，则上式可改写成 （845）式中由式（635）可写成另外由可写出 （846）将式（844）、（845）和（846）合并可写出同步电机的状态方程组，它共含有七个状态变量，方程组的标准形式亦为式中，为非线性项。83 状态方程的数值解法 一、线性定常系统状态方程如果是线性定常系统而非时变系统，通常有拉普拉斯变换求其解析解，此时状态方程的形式为 （847）将式（847）的左右两方进行拉氏变换，可得 （848）式中为单位矩阵，和分别为列

15、向量和的拉氏变换；是的初始值所构成的列向量。式（848）可整理为 （849）于是解答为 （850）式中称为状态转移矩阵，其定义为 （851）由式（850）可见，只要求出系统的状态方程矩阵再经过一些矩阵运算，即可求出。然后再取的拉氏反变换，即可解出向量，即 （852）下面举一实例说明此种解法。【例 81】 已知某一电系统的运动方程为式中为角位移，为外加电流源。设系统的初始条件为零，为单位阶跃函数，试用状态变量法来求解此运动方程。【解】 设，取和为状态变量，则上述运动方程可改写为或于是因 ，又为单位阶跃函数，故于是即将展开成部分分式，可得故 （853）同理将展开成部分分式，可得变换后可得 （854

16、）也可不必将展开成部分分式再进行反变换而直接对式（853）的求其对时间的一阶导数即可得出式（854）式。二、非线性系统和线性时变系统上述的线性定常系统，较易求出其解析解。通常遇到的多是非线性系统和线性时变系统的状态方程 （855）要求出其解析解比较困难，有时是不可能的。从数学观点看，含有n个装态变量的状态方程（855），相当于在n维坐标空间中，由n个状态变量组成的一个矢量。随着时间的推移，此矢量随之变化，而其末端将描绘出一条轨迹。这条轨迹曲线通常可采用数值积分法求解。所谓数值解法，就是找出微分方程的解在一系列离散节点上的近似值。设任意两个相邻节点间的距离均相同，且为，称为步长，可得或，这里。这

17、种已知的初值条件的微分方程数值解得基本点是：求解过程是“进步式”或“递推式”的，即顺着节点的排列次序一步一步地向前推进。例如最简单的解一阶常微方程的欧拉法，其计算格式是，它的含义如图83所示。83欧拉法示意图由图可见，在时的变量值求得后（设精确值称为）应用给定的方程算出在和值下对应的导数值，就可以确定下一个瞬间时的变量值了；换言之，数值积分法是以一组短直线组成的折线来近似代替所研究的曲线，所以它的计算结果是近似的。计算的精度取决于所选用的计算方法和步长的大小。在电机工程中，为了达到足够的精度，最常用的数值积分法有两种：四阶 龙格库塔（Runga-Kutta）法和预估校正法。四阶龙格库塔法的要点

18、是：在每一积分分布内计算四次斜率，再算出加权平均的斜率，然后确定下一瞬间的变量值。其计算公式为 (856)其中 （857）龙格库塔法的主要优点是计算时只要用到值，它是一种一步法，在计算中要改变步长比较方便，而且它具有较高的精度四阶龙格库塔法的截断误差为。但是，在每一时间步内需要计算四次斜率，故计算时间较多。如果改取放大的步长，则其精度即迅速降低。此外，这种方法的稳定性较差，计算机运算时如步长过大会数值溢出，即出现不稳定现象，因此应用范围受到限制。预估校正法则是一种多步法，它利用前面几步中算得的变量值及其一次导数值，在时预估计时的变量值，以后再用校正公式来改进预估值的精度，得出。这种方法的精度取

19、决于每步所使用的前面的和值数目的多少，以及这些值中每一个值带来的相对影响。常用的阿达姆班希福斯（Adams-Bashforth）预估校正公式为 （858）式中为在时的预估值，为其校正值。预估校正法的主要优点在于有较高的精度和稳定性，而且每步只需计算一次斜率，其计算工作量比四阶龙格库塔法少，而截断误差同阶。但是每一步的预估值与前面好几步的数值有关，计算值时要用到、和的值。当所积函数在某段发生突变时，预估值就会有相当大的误差，需要采取另外的弥补措施。84 计算实例【例82】 有一台他励直流电动机，其参数为：，。外施电压为，负载转矩为。试用状态变量法和计算机对此电动机启动过程中的，和求解。【解】 式

20、（816）为他励直流电动机的状态方程，重抄如下：下面用龙格库塔法进行数值积分求解。计算机程序框图如图84所示。设取开关合闸瞬间为，则初始条件为：，。取计算步长为。用FORTRAN语言编制源程序如下。在该程序中，对应的DY(1)DY(4)是这四个变量对时间的导数值。主程序DIMENSION Y (4), DY (4), B (4), C (4)Y (1) =0.0Y (2) =0.0Y (3) =0.0Y (4) =0.0H=0.1CALL RK1 (4, H, 0, Y, DY, B, C)WRITE(5，10)Y, DY10 FORMAT (1X, 4F 15.7)11 CALL RK1 (

21、4, H, 1, Y, DY, B, C)WRITE (5, 10) YIF (Y (1), LT, 15.0) GOTO 11STOPEND子程序 1SUBROUTINE RK1 (N, H, L, Y, DY, B, C)DIMENSION A (4), Y (N), DY (N), B (N), C (N)A (1) = H/2.0A (2) = H/2.0A (3) = HA (4) = HIF (L) 2, 2, 42 DO 3 I=1 ,N3 B (I) = Y (I)CALL F (N, Y, DY)RETURN4 DO 5 K=1, 3DO 6 I = 1, NC (I) =

22、B (I) + A (K)*DY (I)6 Y (I) = Y (I) + A (K+1)*DY (I)/3.05 CALL F (N, C, DY)DO 7 I=1, N7 Y (I) = Y (I) + H*DY (1)/6.0GO TO 2END子程序 2SUBROUTINE F (N, Y, DY)DIMENSION Y (N), DY (N)DY (1) =1.0DY (2) = (220.0-0.1095*Y (4)*Y (3)-0.316*Y (2)/0.142DY (3) = (220.0-178.0*Y (3)/42.72DY (4) = (0.1095*Y (3)*Y (2

23、)-0.00184*Y (4)-4.53)/0.42RETURNEND运用上述程序经电子计算机计算后得到的电枢电流，磁场电流和转速，绘制出曲线如图85所示。最后应当指出，图85中的和的形状虽是典型的但并非唯一的。当点击的参数变动较大时，电流和转速也可能出现振荡。例如当，就会出现这种情况，直到启动15秒末，振荡虽已衰减很多，电流仍在20±160A而转速在800±120rand/s作振荡。图85 他励直流电动机在负载下启动的、和曲线图84 他励直流电动机启动计算用的程序图【例 83】一台三相四极感应电动机，2.2kW，380V，4.81A，1430r/min，定子Y接法。定转子

24、绕组的阻抗参数为：，。感应电动机驱动一他励直流发电机供给纯电阻负载。机械负载制动转矩假定与转速成正比，,为机组的角速度，其单位为rad/s。机组的转动惯量。旋转阻力系数为。要求用综合矢量表示的状态方程来计算电机接上三相对称380V电压从静止状态启动时的动态特性，和。假定，在启动过程中参数都不变。【解】 感应电动机用综合矢量表示的状态方程见式（840）。由，可算得，。再加上给定的数据，可求出状态方程的各项系数。由于FORTRAN语言和四阶龙格库塔法也适用于所要求的复数运算，可采用它们编制程序。考虑到求解的收敛性和稳定性，计算中取步长。程序中，。对应的DY(1)DY(4)是这四个变量对时间的导数。

25、YT为电磁转矩，YR为转速。主程序COMPLEX Y (4), DY (4), B (4), C (4), TY (1) = (0.0, 0.0)Y (2) = (0.0, 0.0)Y (3) = (0.0, 0.0)Y (4) = (0.0, 0.0)H = 0.001CALL CRK1 (4, H, 0, Y, DY, B, C)WRITE (5, 10) Y , DY10 FORMAT (1X, 4F12.4)20 CALL CRK1 (4, H, 1, Y, DY, B, C)YR = 15*Y (4)/3.14159YA = REAL (Y (2)CALL CTORQU (T, Y(

26、2), Y (3)YT = REAL (T)X = REAL (Y (1)WRITE (5, 25) X, YA, YR, YT25FORMAT (1X.T=,F9.6,I=,F9.5,M=,F10.4)IF (REAL (Y(1).LT.0.3) GOTO 20PAUSEEND子程序1SUBROUTINE CRK1 (N, H, L, Y, DY, B, C)COMPLEX Y (N), DY (N), B (N), C (N)DIMENSION A (4)A (1) = H/2.0A (2) = H/2.0A (3) = HA (4) = HIF (L) 30, 30, 50 30DO 4

27、0 I = 1, N40B (I) = Y (L)CALL CF (N, Y, DY)RETURN50DO 60 K = 1,3DO 70 I = 1, NC (I) = B (I) +A (K)*DY (I)70Y (I) = Y (I) +A (K+1)*DY (I)/3.0CALL CF (N, C, DY)60CONTINUEDO 80 I = 1, N80Y (I) = Y (I)+H*DY (I)/6.0GOTO 30END子程序 2SUBROUTINE CF (N, Y, DY)COMPLEX Y (N), DY (N)DY (1) = (1.0, 0.0)DY (2) = -1

28、27.0*Y (2) +81.4*Y (3)-9.53*(0.0, 1.0)*Y (4)*Y (2) -10.08*(0.0, 1.0)*Y (4)* Y (3) +11790* c o s (314*REAL (Y (1) - 0.4712 ) + 11790*(0.0, 1.0)*SIN (314*REAL (Y (1) 0.4712)DY (3) = 120.1 * Y (2) 85.08 *Y (3) +9.96 *( 0.0 , 1.0)*Y (4) *Y (2) + 10.53 * (0.0 , 1.0 ) *Y (4) * Y (3) -11160 *cos (314*REAL (Y (1)-0.4712)-11160*(0.0,1.0)*SIN (314*REAL (