Matlab符号代数运用——电机静态稳定性分析

模型前提假设

  （1）电机铁磁部分的磁路为线性，不计饱和、剩磁、磁滞和涡流的影响；
  （2）电机的定转子齿槽效应用气隙的卡式系数表征，即电机的定转子表面光滑；
  （3）定子三相绕组对称，转子结构分别对直轴和交轴对称。
  （4）不考虑阻尼绕组且发电机无励磁调节器。

数学模型

  通过公式的推到的简化，可以得到如下的特征方程（具体的推导过程可以参考《电机瞬态分析一书》）。
$$D(p)=\left[\begin{array}{ccc}{Z}_d(p)& -x_q(p)& -(u_{q0}+{\psi }_{q0}p)\\ x_d(p)& Z_q(p)& u_{d0}+{\psi }_{d0}p\\ i_{q0}{x}_d(p)-{\psi }_{q0}& {\psi }_{d0}-i_{d0}{x}_q(p)& -H{p}^2\end{array}\right]$$
  当这个特征方程的根均为负值实根时，各微变量将按相应的时间指数函数一直衰减至零，即电机转子将不经任何摇摆而一直回到原来的稳态位置；若特征方程的根均为实根，且出现正值实根时，上述微变量的指数项中将出现有正的时间指数，使其随时间一直增大，从而转子将不经过任何摇摆而偏离原稳态位置，直至脱离同步；如果特征方程的根中只有复数根，但不存在任何正值实部复根时，各微变量将为周期性变化其幅值随时间按指数衰减，此时电机转子将经过振荡逐步回到原位置，如果其复数的实部为正值时，各微变量仍为周期性变化，但其幅值将随时间指数增加，转子经过振荡幅值不断增大而脱离同步。
  对于不考虑阻尼绕组的同步电机来说，且此时满足下述公式：
$$\{\begin{array}{l}{Z}_d(p)=\frac{x_d+x_d\prime T_{d0}p}{1+T_{d0}p}\\ Z_q(p)=x_q\end{array}$$
  此时$D(P)$是一个含有p的五次方的方程，如下所示：
$$D(p)=a_5{p}^5+a_4{p}^4+a_3{p}^3+a_2{p}^2+a_1{p}^1+a_0$$
  可以根据赫尔维兹判据可知，电机的稳定条件为$\Delta$为正定矩阵，其中$\Delta$的表达式如下所示：
$$△=\left[\begin{array}{ccccc}{a}_4& a_2& a_0& 0& 0\\ a_5& a_3& a_1& 0& 0\\ 0& a_4& a_2& a_0& 0\\ 0& a_5& a_3& a_1& 0\\ 0& 0& a_4& a_2& a_0\end{array}\right]$$
  在实际仿真的过程中，由于计算过程过于复杂，故而考率借助$Matlab$中的符号计算得到需要的表达式。

代码思路

  根据上述表达式首先需要求解$D(p)$，但是应当注意到$x_d(p)$的分母问题，故而实际的$Matlab$的代码应该如下所示：

D=[xd*p+Rs,-xq,-(uq0+Wq0*p);
    xd,xq*p+Rs,ud0+Wd0*p;
    iq0*xd-Wq0,Wd0-id0*xq,-p*p*H]
collect((1+p*13.1*2*50*pi)*det(D),p)

  即先求取$D(p)$的行列式结果后，乘以$Z_d(p)$的分母，之后再用$collect$函数以$p$的阶次从高到低排列。之后按照赫尔维兹判据，如果稳定就打点，不稳定就空白可以画出$E_m=1$和$E_m=2$时电机可以静态稳定的图像如下所示：


  为了追求图片的美观，可以通过赫尔维兹判据，将可以稳定的点赋值为$1$，将不可稳定的点赋值为$0$，之后通过等高线绘图得到最终的结果如下所示：