向前欧拉公式 matlab_系统的数学模型与matlab计算

一. 基本知识总结

1.线性系统满足叠加性和齐次性；

2.对呈现非线性特性的元件，通常使用“小信号”方法将这些元件线性化，这也是对电子线路和晶体管进行线性化等效处理的惯用方法；（在工作点处泰勒级数展开 ）

3.通过拉普拉斯变换 ，可以得到特征方程、零点和极点。零点和极点都是特殊的频率点，在极点处Y(s)为无穷大；在零点处Y(s)为零。可以用图示法来表示零点和极点在复频域s平面上的分布，零-极点分布图刻画了系统时间相应的瞬态特性。

4.在对拉普拉斯变换进行处理的过程中，可以使用留数定理。

5.终值定理的使用条件是：Y(s)不能在虚轴和右半平面上存在极点，也不能在原点处存在多重极点，则可以使用公式：

6. 二阶系统的常用表达式与基本参数分析。对于一个二阶系统（质量块-弹簧-阻尼器系统）可以改写为

, 其中，

为无量纲的

阻尼系数，

为系统的

固有（自然）频率。

通过求解特征方程的根，

，则可以分析：

(1)当

>1时，有两个实根，系统为

过阻尼系统；

(2)当

<1时，有一对共轭复根，系统为

欠阻尼系统；

(3)当

=1时，有两个相等的复实根，系统为

临界阻尼系统；

7. 使用框图模型和信号流图模型来描述系统。

二. 基于matlab的仿真计算

1. 动态方程解析解

对于质量块-弹簧-阻尼器系统，质量块的位移响应y(t)的零输入动态响应为：

这个是由系统的微分方程直接求解得到的解析解：

当y(0)=0.15m;

=

rad/s;

在matlab中输入命令：

%   Unforced Response to an Initial Condition

%   Generates Figure 2.41: Spring-mass-damper unforced response.

%

y0=0.15;

wn=sqrt(2);

zeta=1/(2*sqrt(2));

t=[0:0.1:10];

%

c=(y0/sqrt(1-zeta^2));

y=c*exp(-zeta*wn*t).*sin(wn*sqrt(1-zeta^2)*t+acos(zeta));

%

% Boundary Function

%

bu=c*exp(-zeta*wn*t);bl=-bu;

%

% Plot Unforced Response 

%

plot(t,y,t,bu,'k--',t,bl,'k--'),grid

xlabel('Time (sec)'),ylabel('y(t) (meters)')

legend(['omega_n=',num2str(wn),'    zeta=',num2str(zeta)])

2. 多项式的求解与重组

如求解

的解

使用roots求解根，以及使用poly来重构多项式。

p = [1 3 0 4 ];
r = roots(p)

p = poly(r)

3. 多项式的乘积运算与求解多项式的值

使用conv来展开多项式

当变量给定时，使用polyval来求解多项式的值

%   Using conv and polyval to multiply and evaluate 

%   the polynomials  (3 s^2 + 2 s + 1 ) ( s + 4).

%

p=[3 2 1]; q=[1 4];

n=conv(p,q)

%

value=polyval(n,-5)

4. 求解传递函数/状态空间模型、零极点

使用tf求解传递函数；使用ss求解状态空间模型；

使用pole求解传递函数的极点；使用zero求解传递函数的零点。

%   Using tf to create transfer functions.

%

num1=[10]; den1=[1 2 5];

sys1=tf(num1,den1)

%

num2=[1]; den2=[1 1];

sys2=tf(num2,den2)

%

sys=sys1+sys2


p = pole(sys)

z = zero(sys)

5. 绘制零-极点分布图

使用pzmap来绘制系统的零极点图

%    Transfer function example for G(s) and H(s).

%

numg=[6 0 1]; deng=[1 3 3 1]; sysg=tf(numg,deng);

z=zero(sysg)

%

p=pole(sysg)

%

n1=[1 1]; n2=[1 2]; d1=[1 2*i]; d2=[1 -2*i]; d3=[1 3];

numh=conv(n1,n2); denh=conv(d1,conv(d2,d3));

sysh=tf(numh,denh)

%

sys=sysg/sysh

%

% Generate pole-zero map in Figure 2.47

%

pzmap(sys)

注：实际物理系统中极点个数必然大于或等于零点的个数，上图说明了存在多重极点（或多重零点）

6. 传递函数的串联、并联和反馈

使用series实现传递函数串联；

使用parallel实现传递函数的并联；

使用feedback实现反馈；

[sys]=feedback(sys1,[1],sign):单位反馈，sign默认为-1，表示负反馈，+1表示正反馈；

%   Application of the series function.

%

numg=[1];deng=[500 0 0]; sysg=tf(numg,deng);

numh=[1 1];denh=[1 2]; sysh=tf(numh,denh);

sys =series(sysg,sysh);

sys



%    Application of the feedback function with unity feedback.

%

numg=[1]; deng=[500 0 0]; sys1=tf(numg,deng);

numc=[1 1]; denc=[1 2];  sys2=tf(numc,denc);

sys3=series(sys1,sys2);

sys=feedback(sys3,[1])


%   Application of the feedback function with nonunity feedback.

%

numg=[1]; deng=[500 0 0]; sys1=tf(numg,deng);

numh=[1 1]; denh=[1 2]; sys2=tf(numh,denh);

sys=feedback(sys1,sys2);

sys

7. 重复使用上述的命令则可以实现多回路化简

%   Multi-loop block reduction.

%

ng1=[1]; dg1=[1 10]; sysg1=tf(ng1,dg1);

ng2=[1]; dg2=[1 1]; sysg2=tf(ng2,dg2);

ng3=[1 0 1]; dg3=[1 4 4]; sysg3=tf(ng3,dg3);

ng4=[1 1]; dg4=[1 6]; sysg4=tf(ng4,dg4);

nh1=[1 1]; dh1=[1 2]; sysh1=tf(nh1,dh1);

nh2=[2]; dh2=[1]; sysh2=tf(nh2,dh2);

nh3=[1]; dh3=[1]; sysh3=tf(nh3,dh3);

sys1=sysh2/sysg4;

sys2=series(sysg3,sysg4);

sys3=feedback(sys2,sysh1,+1);

sys4=series(sysg2,sys3);

sys5=feedback(sys4,sys1);

sys6=series(sysg1,sys5);

sys=feedback(sys6,[1])

8. 多重回路化简之后，必须消除公因式之后，得到真正的传递函数

使用 minreal 完成零-极点对消

%   Application of the minreal function.

%

num=[1 4 6 6 5 2];

den=[12 205 1066 2517 3128 2196 712];

sys1=tf(num,den);

sys=minreal(sys1)

9. 电力牵引电机控制系统的传递函数与单位阶跃响应

传递函数使用上述的tf、series、feedback可以得到；

单位阶跃响应使用 step可以得到。

%   Application of the minreal function.

%

num=[1 4 6 6 5 2];

den=[12 205 1066 2517 3128 2196 712];

sys1=tf(num,den);

sys=minreal(sys1)



%   Step response of the traction motor wheel velocity.

%

num=[5400]; den=[2 2.5 5402];

sys=tf(num,den);

t=[0:0.005:3];

[y,t]=step(sys,t);

plot(t,y),grid

xlabel('Time (sec)')

ylabel('Wheel velocity')