常用的simulink控制工具箱的函数

常用的simulink控制工具箱的函数

1.构建基本的模型函数1.模型的表示(1)传递函数分子/分母多项式模型(2)传递函数零极点增益模型(3)状态空间模型：2.系统建模(1)串联(2)并联(3)反馈 2.时域响应1.脉冲响应：2.阶跃响应：3.一般响应：4.系统的瞬态性能指标 3.频率特性的分析(1)画$Nyquist$图(2)$Bode$图的绘制(3)频率特征量 4.系统稳定性分析(1)系统特征根(2)系统的相对稳定性 5.系统校正(1)相位超前校正(2)相位滞后校正(3)相位滞后超前校正(4)$PID$校正(5)$Matlab$代码 6.相轨迹图的绘制(1)四阶$Runge-Kutta$法(2) 相轨迹图的绘制 7.线性离散系统(1)$Z$变换的定义(2)$Z$变换的性质1.线性性质：2.延迟定理：3.超前定理：4.初值定理：5.终值定理： (3)$Z$变换表(4)$Matlab$函数(5)$Matlab$实现代码

1.构建基本的模型函数

1.模型的表示

​ 模型的表示有3种基本形式：

传递函数分子/分母多项式模型传递函数零极点增益模型状态空间模型

(1)传递函数分子/分母多项式模型

$$G(s)=\frac{b_m{s}^m+b_{m-1}{s}^{m-1}+...+b_{1}s+b_0}{a_n{s}^n+a_{n-1}{s}^{n-1}+...a_{1}s+a_0}$$

matlab代码： $$\begin{gathered} num=[b_m,b_{m-1},...b_1,b_0],den=[a_n,a_{n-1},...a_1,a_0] \\ G=tf(num,den); \end{gathered}$$

(2)传递函数零极点增益模型

$$G(s)=K\frac{(s-z_0)(s-z_1)...(s-z_m)}{(s-p_0)(s-p_1)...(s-p_n)}$$ matlab代码： $$\begin{gathered} z=[z_0,z_1,...z_m] \\ p=[p_0,p_1,...p_n] \\ k=[K] \\ G=zpk(z,p,k) \end{gathered}$$

复杂函数用$conv(x1,x2)$实现连个向量卷积，用来求多项式乘法。

(3)状态空间模型：

$$\begin{gathered} \dot{X}=AX+Bu \\ Y=CX+Du \end{gathered}$$

matlab代码： $$ss(A,B,C,D)$$

2.系统建模

(1)串联

$$\begin{gathered} T(s)=\frac{Y(s)}{R(s)}=\frac{num}{den}=G_1(s)G_2(s) \\ {G}_1(s)=\frac{num1}{den1} \\ {G}_2(s)=\frac{num2}{den2} \\ [num,den]=series(num1,den1,num2,den2); \end{gathered}$$

(2)并联

$$\begin{gathered} G_1(s)=\frac{num1}{den1} \\ {G}_2(s)=\frac{num2}{den2} \\ T(s)=\frac{Y(s)}{R(s)}=\frac{num}{den} \\ [num,den]=parallel(num1,den1,num2,den2) \end{gathered}$$

(3)反馈

$$\begin{gathered} sign:+1:正反馈，-1:负反馈 \\ G(s)=\frac{num1}{den1} \\ H(s)=\frac{num2}{den2} \\ T(s)=\frac{Y(s)}{R(s)}=\frac{num}{den} \\ [num1,den1]=feedback(num1,den1,num2,den2,sign) \end{gathered}$$ 状态空间模型与传递函数的转化：$[den,num]=ss2tf(F,C,G,D)$。将有理多项式转换成空间状态模型：$[F,C,G,D]=tf2ss(num,den);$

2.时域响应

1.脉冲响应：

$y$:输出响应，$T$:仿真时间，$x$:状态响应（状态空间模型），$sys$:模型，$t$:仿真时间 $$[y,T,x]=impulse(sys,t);$$

2.阶跃响应：

$y$:输出响应，$T$:仿真时间，$x$:状态响应（状态空间模型），$sys$:模型，$t$:仿真时间 $$[y,T,x]=step(sys,t);$$

3.一般响应：

$y$:输出响应，$T$:仿真时间，$x$:状态响应（状态空间模型），$sys$:模型，$t$:仿真时间,$u$:输入 $$[y,T,x]=lsim(sys,u,t);$$

4.系统的瞬态性能指标

求$G(s)=\frac{50}{0.05s^2+(1+50\tau )s+50}$在不同的$\tau$时的单位阶跃响应：

上升时间函数:

function trl = trlf( x,yss,t0 ) %trl 计算上身时间 % 此处显示详细说明 r = 1; while x(r)<yss r = r+1; end trl = (r-1)*t0; end

峰值时间函数：

function tpl = tplf( x,t0 ) %tpl 计算峰值时间 % 此处显示详细说明 [ymax,tp] = max(x); tpl = (tp-1)*t0; end

最大超调量和调整时间函数：

function [tsl,mpl] = mplf( x,t0,yss ,dlta) %mplf 计算1001个采样点时的调整时间和最大超调量 % 此处显示详细说明 [ymax,tp] = max(x); mpl = (ymax-yss)/yss; s = 1001; while (x(s)>(1-dlta))&&(x(s)<1+dlta) s = s-1; end tsl = (s-1)*t0; end

主函数：

t = 0:0.001:1; yss = 1; dta = 0.02; nG = [50,50,50]; tao = [0,0.0125,0.025]; trl = zeros(3,1);%上升时间 tpl = zeros(3,1);%峰值时间 mpl = zeros(3,1);%最大超调量 tsl = zeros(3,1);%调整时间 for i = 1:3 dG(i,:) = [0.05,1+50*tao(i),50]; G(i) = tf(nG(i),dG(i,:)); end out1 = step(G(1),t); out2 = step(G(2),t); out3 = step(G(3),t); %以下是计算上升时间 trl(1) = trlf( out1,yss,0.001 ); trl(2) = trlf( out2,yss,0.001 ); trl(3) = trlf( out3,yss,0.001 ); %以下是计算峰值时间 tpl(1) = tplf( out1,0.001 ); tpl(2) = tplf( out2,0.001 ); tpl(3) = tplf( out3,0.001 ); %以下是计算最大超调量和调整时间 [tsl(1),mpl(1)] = mplf(out1,0.001,yss,dta); [tsl(2),mpl(2)] = mplf(out2,0.001,yss,dta); [tsl(3),mpl(3)] = mplf(out3,0.001,yss,dta); %输出 disp('tao = 0 时的上身时间 峰值时间 最大超调量 调整时间： '); disp([trl(1),tpl(1),mpl(1),tsl(1)]); disp('tao = 0.0125 时的上身时间 峰值时间 最大超调量 调整时间： '); disp([trl(2),tpl(2),mpl(2),tsl(2)]); disp('tao = 0.025 时的上身时间 峰值时间 最大超调量 调整时间： '); disp([trl(3),tpl(3),mpl(3),tsl(3)]);

结果：

>> abc tao = 0 时的上身时间 峰值时间 最大超调量 调整时间： 0.0640 0.1050 0.3509 0.3530 tao = 0.0125 时的上身时间 峰值时间 最大超调量 调整时间： 0.0780 0.1160 0.1523 0.2500 tao = 0.025 时的上身时间 峰值时间 最大超调量 调整时间： 0.1070 0.1410 0.0415 0.1880

3.频率特性的分析

(1)画$Nyquist$图

$re:$时频特性，$im:$虚频特性，$w:$频率范围，$sys:$模型，$w$可选频率 $$[re,im,w]=nyquist(sys,w);$$ 例如：绘制$G(s)=\frac{24(0.25s+0.5)}{(5s+2)(0.05s+2)}$的$Nyqusit$图：

k = 24; numG = k*[0.25 0.5]; denG = conv([5,2],[0.05,2]); sys = tf(numG,denG); [m,n] = size(re); [re,im] = nyquist(sys); for i = 1:n re1(i) = re(1,1,i); im1(i) = im(1,1,i); end plot(re1,im1);grid on; xlabel('实部'); ylabel('虚部'); title('系统Nyquist图');

(2)$Bode$图的绘制

$mag:$幅频特性，$phase:$相频特性，$w:$频率范围，$sys:$模型 $$[mag,phase,w]=bode(sys,w);$$ 例如：绘制$G(s)=\frac{24(0.25s+0.5)}{(5s+2)(0.05s+2)}$的$Bode$图：

k = 24; numG = k*[0.25 0.5]; denG = conv([5,2],[0.05,2]); sys = tf(numG,denG); [m,n] = size(re); w = logspace(-2,3,100); bode(sys,w); title('系统Nyquist图');

(3)频率特征量

k = 24; numG = k*[0.25 0.5]; denG = conv([5,2],[0.05,2]); sys = tf(numG,denG); w = logspace(-2,4,100); [Gm,Pm,w]=bode(sys,w); [m,n] = size(Gm); for i = 1:n Gm(i) = Gm(1,1,i); Pm(i) = Pm(1,1,i); end [Mr,k] = max(Gm); disp('协震峰值是：'); disp(20*log10(Mr)); disp('谐振频率是：'); disp(w(k)); n = 1; while 20*log10(Gm(n))>=-3 n = n+1; end disp('截至频率是：'); disp(w(n));

结果：

协震峰值是： 9.5398 谐振频率是： 0.0100 截至频率是： 3.5112

4.系统稳定性分析

(1)系统特征根

已知特征方程$a_n{s}^n+a_{n-1}{s}^{n-1}+...+a_{1}s+a_0$的特征系数是$[a_n,a_{n-1},...a_1,a_0]$，求根的方法： $$roots([a_n,a_{n-1},...a_1,a_0]);$$

(2)系统的相对稳定性

$sys:$系统模型，$Gm:$幅值裕度，$P_m:$相位裕度，$W_{cg}:$相位穿越频率，$W_{cp}:$幅值穿越频率。 $$\begin{gathered} [G_m,P_m,W_{cg},W_{cp}]=margin(sys); \\ [mag,phase,w]=bode(sys);\to [G_m,P_m,W_{cg},W_{cp}]=margin(mag,phase,w) \end{gathered}$$

5.系统校正

(1)相位超前校正

$$\begin{gathered} T=R_1{C}_1 \\ \alpha =\frac{R_1}{R_1+R_2}<1 \\ G(s)=\frac{U_0(s)}{U_i(s)}=\frac{1+Ts}{1+\alpha Ts} \end{gathered}$$

(2)相位滞后校正

$$\begin{gathered} \beta =\frac{R_3+R_4}{R_4}>1 \\ T=R_4{C}_2 \\ G(s)=\frac{1+Ts}{1+\beta Ts} \end{gathered}$$

(3)相位滞后超前校正

$$\begin{gathered} T_1=R_1{C}_1,T_2=R_2{C}_4 \\ \beta =\frac{R_1+R_2}{R_2}>1 \\ G(s)=\frac{(T_{1}s+1)(T_{2}s+1)}{(\frac{T_1}{\beta }s+1)(\beta T_{2}s+1)} \end{gathered}$$

(4)$PID$校正

$$\begin{gathered} T_i=R_1{C}_1+R_2{C}_2 \\ {T}_d=\frac{R_1{C}_1{R}_2{C}_2}{R_1{C}_1+R_2{C}_2} \\ {K}_p=\frac{R_1{C}_1+R_2{C}_2}{R_1{C}_2} \\ G(s)=K_p(1+\frac{1}{T_{i}s}+T_{d}s) \end{gathered}$$

(5)$Matlab$代码

function y = PID_fun(num,den,kp,ki,kd,ts,numbers,yd,M) %PID_fun 阶跃信号 % num:连续函数分子 % den:连续函数分母 % k_p,k_i,k_d:PID参数 % ts:采样时间 % numbers:采样点数目 % yd:输入信号 % M:阈值 % clc; % clear; sys = tf(num,den); dsys = c2d(sys,ts,'z'); [numz,denz] = tfdata(dsys,'v'); u_1 = 0;u_2 = 0;u_3 = 0; y_1 = 0; y_2 = 0; y_3 = 0; x = [0 0 0]'; error_1 = 0; error_2 = 0; for k = 1:1:numbers time(k) = k*ts; du(k) = kp*x(1)+kd*x(2)+ki*x(3); u(k) = u_1 + du(k); if u(k)>=M u(k) = M; end if u(k)<=-M u(k) = -M; end y(k) = (-denz(2)*y_1 - denz(3)*y_2+numz(2)*u_1+numz(3)*u_2)/denz(1);%z变换延迟定理 error = yd(k) - y(k); u_3 = u_2;u_2 = u_1;u_1 = u(k); y_3 = y_2;y_2 = y_1;y_1 = y(k); x(1) = error - error_1; x(2) = error - 2*error_1+error_2; x(3) = error; error_2 = error_1; error_1 = error; end figure; plot(time,yd,'r',time,y,'b','linewidth',2); xlabel('时间(s)');ylabel('y_d,y'); grid on title('增量式PID响应曲线'); legend('理想位置信号','位置追踪'); end

调用

>> y = PID_fun([40 10],[1 50 0],8,0.1,10,0.001,1000,2*ones(1000,1),10);

反馈校正和顺馈校正此处就不再赘述了。

6.相轨迹图的绘制

(1)四阶$Runge-Kutta$法

​ 设$t$为自变量时间，$y$为因变量，微分方程的形式$\dot{y}=f(t,y)$。${}^{\prime }odefu{n}^{\prime }$:包含微分方程(组)的$M$文件，$ts:$自变量取值，$y_0:$变量$y$的初值，$p_i:$文件内的附加参数。 $$[t,y]=ode45{(}^{\prime }odefu{n}^{\prime },t_s,y_0,[],p_1,p_2,...);$$

(2) 相轨迹图的绘制

​ 做出$\ddot{x}+0.5\dot{x}+2x+x^2=0$在初值为$(-2.2,0.5{)}^T$时的相轨迹图： ​ 令$x_1=x,x_2=\dot{x_1}$: $$\begin{gathered} \dot{x_1}=x_2 \\ \dot{x_2}=-0.5x_2-2x_1-x_1^2 \end{gathered}$$ 函数：

function y = odefun1( t,x,p1 ) %UNTITLED 此处显示有关此函数的摘要 % 此处显示详细说明 y = [x(2);-p1*x(2)-2*x(1)-x(1)^2]; end

主函数：

t = 0:0.01:100; x0 = [-2.2;0.5]; p1 = 0.5; [t,y] = ode45(@odefun1,t,x0,[],p1); plot(y(:,1),y(:,2),'linewidth',2); grid on;xlabel('x_1');ylabel('x_2'); title('相轨迹图');

7.线性离散系统

(1)$Z$变换的定义

​ 连续信号$x(t)$,采样输出信号$x^{\ast }(t)$,单位脉冲序列${\delta }_s(t)$,采样周期$T$ $$\begin{gathered} x^{\ast }(t)=\sum _{n=0}^{\infty }x(nT)\delta (t-nT) \\ 当n\ge 0时，对上面进行Laplace变换： \\ L[x^{\ast }(t)]=\sum _{n=0}^{\infty }x(nT){\int }_0^{\infty }\delta (t-nT)e^{-st}dt=\sum _{n=0}^{\infty }x(nT)e^{-snT} \end{gathered}$$ 如果我们让$z=e^{sT}$得到$Z$变换的变换式： $$Z[x(t)]=X(z)=\sum _{n=0}^{\infty }x(nT)z^{-n}$$

(2)$Z$变换的性质

1.线性性质：

$$Z[a{x}_1(t)+b{x}_2(t)]=a{X}_1(z)+b{X}_2(z)$$

2.延迟定理：

$$设Z[x(t)]=X(z),且t<0时，x(t)=0,则:Z[x(t-mT)]=z^{-m}X(z)$$

3.超前定理：

$$设Z[x(t)]=X(z),则：Z[x(t+mT)]=z^m[X(z)-\sum _{k=0}^{m-1}x(kT)z^{-k}]$$

4.初值定理：

$$设Z[x(t)]=X(z),则：x(0)=li{m}_{z\to \infty }X(z)$$

5.终值定理：

$$设Z[x(t)]=X(z),且(z-1)X(z)的全部极点位于单位圆内，则：x(\infty )=li{m}_{z\to 1}[X(z)(z-1)]$$

(3)$Z$变换表

$$\begin{gathered} Z[a^k]=\frac{z}{z-a} \\ Z[a^{k}cosk\pi ]=\frac{z}{z+a} \end{gathered}$$

稳定性分析和校正设计就不描述了。

(4)$Matlab$函数

若$G(z)=\frac{numd(z)}{dend(z)}$,$T$为采样时间，${}^{\prime }zo{h}^{\prime }$为零阶保持，传递函数$G_p(s)=\frac{num(s)}{den(s)}$,且$G(z)$对应着传递函数$G_p(s)$,则： $$[numd,dend]=c2dm(num,den,T{,}^{\prime }zo{h}^{\prime });$$ 而如果求得逆的话： $$[num,den]=d2cm(numd,dend,T{,}^{\prime }zo{h}^{\prime });$$ 如果要求任意输入响应，$y$:输出响应，$x:$状态响应，$u:$输入，$G(s)=\frac{num}{den}$,$k:$指定的采样数。 $$[y,x]=dlsim(num,den,u,k);$$

(5)$Matlab$实现代码

num = 1;den = [1,1,0]; T = 1; [nz,dz] = c2dm(num,den,T,'zoh'); printsys(nz,dz,'z');

输出：

num/den = 0.36788 z + 0.26424 ------------------------ z^2 - 1.3679 z + 0.36788