实例1-电机速度控制(Matlab实现)

Tommy！

已于 2024-07-22 00:30:41 修改

阅读量1.8k

点赞数 42

文章标签：数学建模 matlab

于 2024-07-21 13:37:37 首次发布

本文链接：https://blog.csdn.net/weixin_51297959/article/details/140468029

版权

1.数学建模及分析

这是一个直流电机的等效电路，假设系统的输入就是电压 $v$ ，输出是电动机轴的转速 $\dot{\theta }$ 。通过观察这个模型，我们可以先基于基尔霍夫电压定律，建立电压回路方程。

$V=L\frac{\mathrm{d} i}{\mathrm{d} t}+Ri+e$

上式中的e表示永磁体转子旋转在绕组中产生的反电动势，它与转子的转速 $\dot{\theta }$ 成正比，这里我们引入一个反电动势系数 $K_{e}$ ，可以得到

$e=K_{e}\dot{\theta }$

然后基于牛顿第二定律分析转矩的平衡方程，一般来说，直流电机产生的转矩和电流以及磁场的强度成正比，为了简便题目，我们把磁场的强度看作一个常数，那么就可以得到产生的转矩大小为

$T_{1}=K_{t}i$

$K_{t}$ 为电机转矩系 $b$ 数，除了电枢产生的转矩之外，我们的电动机和负载折合到电动机轴上的粘性摩擦,其转矩与转轴的角速度成正比，假设粘性摩擦系数为 $b$ ,那么粘性摩擦的转矩为

$T_{2}=b\dot{\theta }$

我们知道合外力矩等于转动惯量乘以角加速度，因此可以得到

$T_{1}-T_{2}=J\ddot{\theta }$

在国际单位制中，电机转矩系数等于反电动势系数，即 $K_{t}=K_{e}$ ,因此我们都用 $K$ ，至此我们可以合并上述公式得到关于这个直流电机的数学模型，如下

$V=L\frac{\mathrm{d} i}{\mathrm{d} t}+Ri+K\dot{\theta }$ （1） $Ki-b\dot{\theta }=J\ddot{\theta }$ （2）

我们为了后续的控制器设计，我们需要将上述的数学模型进行拉普拉斯变换，可以得到

$V(s)=LsI(s)+RI(s)+K_{e}\theta(s)$ （3） $KI(s)-bs\theta (s)=Js^{2}\theta (s)$

通过约掉中间量 $i$ ，可以得到传递函数

$H(s)=\frac{\dot{\theta} (s)}{V(s)}=\frac{K}{(Js+b)(Ls+R)+K^{2}}$ （4）

在matlab软件中设置好系统的基本参数，然后通过tf函数构造传递函数

J = 0.01;%转动惯量
b = 0.1;%粘性摩擦系数
K = 0.01;%反电动势系数或电机转矩系数
R = 1;%电阻值
L = 0.5;%电感
s = tf('s');
P_motor = K/((J*s+b)*(L*s+R)+K^2)%传递函数

输出的结果为

通过matlab中的线性系统分析器可以初步分析这个开环传递函数的情况，这分析器挺好用的，单机右键选择plot types可以选择输入是阶跃输入，冲激信号输入，看伯德图，零极点分布图以及奈奎斯特图都可以。

linearSystemAnalyzer('step', P_motor, 0:0.1:5);%输入单位阶跃响应

运行代码，得到响应结果

从图中我们可以看到，响应得到的最大值为0.1，并且需要2.07秒的时间达到稳定状态（要观察这些点的方式只需要右键单击图片，然后选择characteristics，再选择setting time,rise time和steady state）。

对于我们的电机速度系统，我们的期望是其稳态误差要小于1%，并且设定时间要小于2秒，超调量小于5%，那么这个开环系统的阶跃响应很明显是不符合要求的。而后续的控制器的目标就是要显著的减少稳态误差（从0.1提高到1左右），略微减少设定时间（从2.07减少到低于2s就可以了），同时做到超调量不超过5%。

2.状态空间法设计控制器

要利用状态空间法设计控制器首先得写出状态空间表达式，在第一章的基础上，我们选择转子的角速度 $\dot{\theta }$ 和电流 $i$ 为状态变量，通过式（1），式（2），我们可以得到状态方程为

$\frac{\mathrm{d} }{\mathrm{d} t}\begin{bmatrix} \dot{\theta} \\ i \end{bmatrix}=\begin{bmatrix} -\frac{b}{J} &\frac{K}{J} \\ -\frac{K}{L} & -\frac{R}{L} \end{bmatrix}\begin{bmatrix} \dot{\theta} \\ i\end{bmatrix}+\begin{bmatrix} 0\\ \frac{1}{L}\end{bmatrix}V$

输出方程为

$y=\begin{bmatrix}1 & 0\end{bmatrix}\begin{bmatrix} \dot{\theta }\\ i \end{bmatrix}$

将代码写成状态方程形式，如下所示。

J = 0.01;
b = 0.1;
K = 0.01;
R = 1;
L = 0.5;
A = [-b/J   K/J
    -K/L   -R/L];
B = [0
    1/L];
C = [1   0];
D = 0;
sys = ss(A,B,C,D)

2.1全状态反馈控制器

将 $x^{T}=\begin{bmatrix} \dot{\theta } & i \end{bmatrix}$ ， $A=\begin{bmatrix} -\frac{b}{J} &\frac{K}{J} \\ -\frac{K}{L} & -\frac{R}{L} \end{bmatrix}$ ， $B=\begin{bmatrix} 0\\ \frac{1}{L}\end{bmatrix}$ ， $C=\begin{bmatrix} 1\\ 0\end{bmatrix}$ ， $D=0$ ,可以把状态空间方程画成如下图，全状态反馈的本质就是在状态变量 $x$ 处进行比例反馈，接下来就是确定 $K_{c}$ 的值。

在确定 $K_{c}$ 的值之前，还需要做的工作是判断这个系统是否具有可控性，而判断的依据就是矩阵 $\begin{bmatrix} B &AB &A^{2} B&... & A^{n-1}B \end{bmatrix}$ 是否为n。

在matlab中应用ctrb函数可以直接得到系统的可控性矩阵

sys_rank = rank(ctrb(A,B))%可控性矩阵的秩

输出结果为2。

再看系统的阶数

sys_order = order(sys)%获取模型阶数

输出结果也是2，因此满足可控性条件。

接下来就是考虑闭环极点的位置放置了，我们知道一个闭环反馈系统的稳定性主要是根据 $\left ( A-BK \right )$ 的特征值位置所决定的，特征值等于闭环系统极点的位置，因此我们需要先决定闭环系统极点的位置。

先假设极点位置为 $-5+i$ , $-5-i$ 。(这个极点位置是我们希望的系统极点位置，然后它的大小是要根据我们的设计要求指标，即我们的设定时间，稳态误差和超调量所决定的，具体计算细节可以参考下图中的步骤。)

假设好希望极点的位置后，我们需要计算 $K_{c}$ 值，理论的计算方式是我们需要求解出状态反馈系统的特征方程然后对应上我们希望极点位置下的特征方程，将其一一对应，然后求出 $K_{c}$ 矩阵中的两个值，这样的话比较麻烦。

在matlab中有一个place函数，可以输入两个希望极点的位置和状态方程的A,B矩阵，然后直接得出反馈增益 $K_{c}$ 的值。

Kc = place(A,B,[p1 p2]);确认反馈增益大小

输出结果为

这样的话，我们就可以得到我们的闭环系统状态空间表示

$\dot{x}=(A-B*Kc)x+Br$

$y=Cx$

上式中的由来，通过5.1节的第一张图就能明白了。

然后就需要去得到闭环系统的阶跃响应，观察图像。

sys_cl = ss(A-B*Kc,B,C,D);%闭环系统
linearSystemAnalyzer('step', sys_cl, 0:0.01:5);%线性系统分析器

输出为

观察图像可以知道，设定时间为1.1s，但是终值只有0.0769，也就是说稳态误差远大于预期，超调量接近于0。因此接下来就通过增加一个前置补偿器减少系统的稳态误差。

2.2增加一个前置补偿器

前置补偿器说白了就是通过放大或者是缩小我们的输入值(成一定倍数)，从而放大我们的输出，进而减少稳态误差。而补偿的数值大小，需要我们根据当前的输出来进行修正。这里给出一个计算补偿值Nbar的函数，函数内容的含义可以看看每句话的解释意思。

function[Nbar]=rscale(a,b,c,d,k)
error(nargchk(2,5,nargin));
%判断输入是否符合要求
nargin1 = nargin;
if (nargin1==2), % 系统形式
[A,B,C,D] = ssdata(a);
K=b;
elseif (nargin1==5), % A,B,C,D matrices
A=a; B=b; C=c; D=d; K=k;
else error('Input must be of the form (sys,K) or (A,B,C,D,K)')
end;
% 计算 Nbar
s = size(A,1);%获取矩阵 A 的行数，假设为状态空间维度
Z = [zeros([1,s]) 1];%创建一个行向量 Z，前 s 个元素为零，最后一个元素为 1
N = inv([A,B;C,D])*Z';%计算增广矩阵 [A,B;C,D] 的逆乘以向量 Z，得到向量 N
Nx = N(1:s);%提取向量 N 的前 s 个元素，赋值给 Nx
Nu = N(1+s);%提取向量 N 的第 s+1 个元素，赋值给 Nu
Nbar=Nu + K*Nx;%计算 Nbar，即 Nu + K * Nx，这是用于控制系统设计中的输入修正值

然后就是具体的计算Nbar的值以及将它加入到闭环系统的阶跃响应当中去。

sys_cl = ss(A-B*Kc,B,C,D);%闭环系统
Nbar = rscale(sys,Kc);%计算补偿值
linearSystemAnalyzer('step', sys_cl*Nbar, 0:0.01:5);

输出的结果为

可以发现通过增加补偿后，系统的稳态误差大大减少了，也没有超调量，设定时间为1.1s，符合我们最一开始的系统设计要求，因此就算成功的通过状态空间法确定了此系统的控制器设计。

3.结束

这次是状态空间表示法比较简单的应用，要分析一切复杂的系统肯定是需要先从基础开始，然后再一步一步提升的，所以后面会进一步选择一些更加复杂的系统来进行控制器的设计。这个简单系统的控制器可以选择PID，根轨迹，频域法等等方法都可以实现，目前我学习的侧重点在状态空间，因此我就只选择了通过状态空间法来实现这个系统的控制器设计。