建立三相异步电机在两相旋转dq坐标系 下的定子数学模型，其中包含电压方程、磁链方程、电磁转矩方程、机械运动方程和机械负载等部分，此电机模型用于后续的机电耦合建模，采用matlab中的ode45编写程序，注意先写主程序再写调用函数将两个部分分开写完，电机参数请使用地铁电机参数

三相异步电机在两相旋转dq坐标系下的定子数学模型可以表示为：

1. 电压方程

$$v_d=R_si_d+L_s\frac{di_d}{dt}-\omega_s L_s i_q$$

$$v_q=R_si_q+L_s\frac{di_q}{dt}+\omega_s L_s i_d$$

其中，$v_d$和$v_q$分别为d轴和q轴上的电压，$R_s$为定子电阻，$L_s$为定子电感，$i_d$和$i_q$分别为d轴和q轴上的电流，$\omega_s$为同步速度。

2. 磁链方程

$$\frac{d\psi_d}{dt}=v_d-\omega_s\psi_q$$

$$\frac{d\psi_q}{dt}=v_q+\omega_s\psi_d$$

其中，$\psi_d$和$\psi_q$分别为d轴和q轴上的磁链。

3. 电磁转矩方程

$$T_e=\frac{3}{2}P\frac{p_s}{p_r}\frac{L_s}{\omega_s}(i_d\psi_q-i_q\psi_d)$$

其中，$T_e$为电磁转矩，$P$为功率因数，$p_s$和$p_r$分别为定子极数和转子极数。

4. 机械运动方程

$$J\frac{d\omega_m}{dt}=T_e-T_l$$

其中，$J$为转动惯量，$\omega_m$为机械角速度，$T_l$为机械负载。

5. 机械负载

$$T_l=B\omega_m+C\omega_m^2+T_{l0}$$

其中，$B$为摩擦系数，$C$为阻尼系数，$T_{l0}$为静摩擦力矩。

根据以上方程，可以在matlab中编写程序，使用ode45求解。程序分为主程序和调用函数两部分，具体代码如下：

主程序：

clc;clear all;close all;
% 定义电机参数
Rs = 1.82;                  % 定子电阻
Ls = 0.0271;                % 定子电感
P = 3;                      % 功率因数
ps = 3;                     % 定子极数
pr = 2;                     % 转子极数
J = 0.02;                   % 转动惯量
B = 0.1;                    % 摩擦系数
C = 0.01;                   % 阻尼系数
Tl0 = 0;                    % 静摩擦力矩
% 定义初始状态
id0 = 0;                    % 初始d轴电流
iq0 = 0;                    % 初始q轴电流
w0 = 0;                     % 初始机械角速度
% 定义仿真时间和步长
tspan = [0 0.2];            % 仿真时间范围
h = 0.0001;                 % 步长
% 调用ode45求解函数
[t,y] = ode45(@(t,y) motor(t,y,Rs,Ls,P,ps,pr,J,B,C,Tl0),tspan,[id0,iq0,w0]);
% 绘图
figure(1)
subplot(3,1,1)
plot(t,y(:,1),'LineWidth',1.5)
xlabel('Time(s)')
ylabel('i_d(A)')
grid on
subplot(3,1,2)
plot(t,y(:,2),'LineWidth',1.5)
xlabel('Time(s)')
ylabel('i_q(A)')
grid on
subplot(3,1,3)
plot(t,y(:,3),'LineWidth',1.5)
xlabel('Time(s)')
ylabel('w(rad/s)')
grid on

调用函数：

function dydt = motor(t,y,Rs,Ls,P,ps,pr,J,B,C,Tl0)
% 电压方程
vd = Rs*y(1) + Ls*(diff(y(1:2))/diff(t)) - y(3)*Ls*y(2);
vq = Rs*y(2) + Ls*(diff(y(2:3))/diff(t)) + y(3)*Ls*y(1);
% 磁链方程
dpsidt(1) = vd - y(3)*y(2);
dpsidt(2) = vq + y(3)*y(1);
% 电磁转矩方程
Te = 1.5*P*ps/pr*Ls/y(3)*(y(1)*dpsidt(2)-y(2)*dpsidt(1));
% 机械运动方程
dwdt = (Te-Tl0-B*y(3)-C*y(3)^2)/J;
% 输出
dydt = [dpsidt(1);dpsidt(2);dwdt];
end

需要注意的是，在调用函数中，磁链方程的导数需要分别计算，因此用一个数组`dpsidt`来存储。另外，在主程序中，使用`subplot`函数将三个变量随时间变化的曲线绘制在同一张图中。最后，运行程序即可得到结果。