西安石油大学电子工程学院：《自动控制理论 Modern Control System》精品课程教学资源（课程设计）指导书

《自动控制理论》课程设计指导书 霍爱清薛朝妹 自动化教研室

《自动控制理论》课程设计指导书 霍爱清 薛朝妹 自动化教研室

目录 1控制系统的数学描述… 1.1微分方程… 1.2传递函数… 1.3状态空间描述… 14模型转换… 2控制系统的校正… …15 2.1单变量系统的两种主要校正方式… 2.2PI、PD、PID校正… 15 2.3串联校正举例…………………… …18 3 MATLAB在自动控制系统中的应用… 3.1概述………… 24 3.2控制系统函数全集 …30 3.3应用实例…… 4 SIMULINK简介 4.1引导 42 SIMULINK模型的构建……… 5自动控制系统设计任务……… 51任务一………… 52任务 53任务 54任务四……………… …54 6附录… 56 6.1时域分析例题及程序………………… 62根轨迹分析例题及程序 6.3频域分析例题及程序 6

目 录 1 控制系统的数学描述…………………………………………………………………………1 1.1 微分方程………………………………………………………………………………1 1.2 传递函数………………………………………………………………………………6 1.3 状态空间描述…………………………………………………………………………9 1.4 模型转换………………………………………………………………………………11 2 控制系统的校正………………………………………………………………………………15 2.1 单变量系统的两种主要校正方式……………………………………………………15 2.2 PI、PD、PID 校正……………………………………………………………………15 2.3 串联校正举例…………………………………………………………………………18 3 MATLAB 在自动控制系统中的应用………………………………………………………24 3.1 概述……………………………………………………………………………………24 3.2 控制系统函数全集……………………………………………………………………30 3.3 应用实例………………………………………………………………………………31 4 SIMULINK 简介………………………………………………………………………………41 4.1 引导……………………………………………………………………………………41 4.2 SIMULINK 模型的构建………………………………………………………………45 5 自动控制系统设计任务………………………………………………………………………50 5.1 任务一…………………………………………………………………………………50 5.2 任务二…………………………………………………………………………………52 5.3 任务三…………………………………………………………………………………53 5.4 任务四…………………………………………………………………………………54 6 附录…………………………………………………………………………………………56 6.1 时域分析例题及程序…………………………………………………………………56 6.2 根轨迹分析例题及程序………………………………………………………………64 6.3 频域分析例题及程序…………………………………………………………………66

第一章控制系统的数学描述 1.1微分方程 111物理系统的微分方程 利用机械学、电学、流体力学和热力学等的物理规律,我们可以得到物理系统的动态 方程。它们通常用常系数线性微分方程来描述 1.1.2数值解 通过拉普拉斯变换和反变换,可得到线性时不变方程的解析解,也可用状态转移矩阵 φ(t)求解。这些分析方法通常只限于常系数的线性微分方程。解析解是精确的,然而通常寻 找解析解是困难的,甚至是不可能的。而数值分析方法直接在时域里求解微分方程,不仅适 用于线性时不变方程,也适用于非线性以及时变微分方程。 MATLAB提供了两个求微分方程数值解的函数,它们采用龙格-库塔( Ode23和ode45分别表示采用2阶和4阶龙格-库塔公式,后者具有更高的精度 n阶微分方程必须化为n个首1的一阶微分方程组,且放入M文件中,以便返回方程 状态变量的导数,下面的例子介绍这些函数的用法。 例1.1对图1-1的机械系统,已知三个量一一拉力、摩擦力、以及弹簧力都影响质量 M的加速度 利用牛顿运动定理,建立系统的力平衡方程式 d x+b-+Kx=f(t) B d. 令x1=X,X2=,有 (t) x(t) dx dtx2 图1-1 =[f(t)-Bx2-Kx1] 设质量M=lkg,摩擦系数B=5Nm/sec,弹簧常数K=25N/m。在=0时刻,施加25N的 拉力。上述方程及已知量在M文件 mechsys m中定义如下 function dot=mechsys(t, x) F=25

1 第一章 控制系统的数学描述 1.1 微分方程 1.1.1 物理系统的微分方程 利用机械学 、电学、流体力学和热力学等的物理规律，我们可以得到物理系统的动态 方程。它们通常用常系数线性微分方程来描述。 1.1.2 数值解 通过拉普拉斯变换和反变换，可得到线性时不变方程的解析解，也可用状态转移矩阵 φ(t)求解。这些分析方法通常只限于常系数的线性微分方程。解析解是精确的，然而通常寻 找解析解是困难的，甚至是不可能的。而数值分析方法直接在时域里求解微分方程，不仅适 用于线性时不变方程，也适用于非线性以及时变微分方程。 MATLAB 提供了两个求微分方程数值解的函数，它们采用龙格-库塔(Runge-kutta)法。 Ode23 和 ode45 分别表示采用 2 阶和 4 阶龙格-库塔公式，后者具有更高的精度。 n 阶微分方程必须化为 n 个首 1 的一阶微分方程组，且放入 M-文件中，以便返回方程 状态变量的导数，下面的例子介绍这些函数的用法。 例 1.1 对图 1-1 的机械系统，已知三个量——拉力、摩擦力、以及弹簧力都影响质量 M 的加速度。 利用牛顿运动定理，建立系统的力平衡方程式 Kx f(t) dt dx B d t d x M 2 2    令 dt dx x1  x, x2  ，有 2 1 x dt dx  设质量 M=1kg，摩擦系数 B=5N/m/sec，弹簧常数 K=25N/m。在 t=0 时刻，施加 25N 的 拉力。上述方程及已知量在 M-文件 mechsys.m 中定义如下： function xdot=mechsys(t, x); F=25; [f(t) Bx Kx ] M 1 dt dx 2 1 2   

M=1B=5K=25 xdot=[x(2);M*(F-B*x(2)-K*x(1) 下面的M文件使用ode23对系统在零初始条件下进行仿真: t0=0: tfinal=3 %时间间隔0~3秒 x0=[0,0]; %零初始条件 tol=0.001; %精度 %如果非零,则打印出每一步的计算值 It, x]=ode23(mechsys, to, tfinal, xO, tol,trace) ubplot(21 1),plot(t, x) title (Time response of mechanical translational system) text(2, 1.2, displacement) text(2,. 2, veloclty') subplot(2 12),plot(d, v); title (velocity versus displacement) xlabel (displacement) ylabel (velocity) subplot(lll) 仿真结果如图1-2 Time response of mechanical translational system 2 0.6 displacement 图1-2

2 M=1;B=5;K=25; xdot=[x(2);1/M*(F-B*x(2)-K*x(1))]; 下面的 M-文件使用 ode23 对系统在零初始条件下进行仿真： t0=0; tfinal=3; ％时间间隔 0～3 秒 x0=[0,0]; ％零初始条件 tol=0.001; ％精度 trace=0; ％如果非零，则打印出每一步的计算值 [t, x]=ode23(’mechsys’,t0,tfinal,x0,tol,trace) subplot(211),plot(t, x); title (’Time response of mechanical translational system’) xlabel (’Time-sec’) text (2,1.2,’displacement’) text (2,.2,’veloclty’) d=x(:,1);v=x(:,2); subplot(212),plot(d,v); title (’velocity versus displacement’) xlabel (’displacement’) ylabel (’velocity’) subplot(111) 仿真结果如图 1-2。 0 0.5 1 1.5 2 2.5 3 -1 0 1 2 3 Time response of mechanical translational system Time-sec displacement veloclty 0 0.2 0.4 0.6 0.8 1 1.2 1.4 -1 0 1 2 3 velocity versus displacement displacement velocity 图 1-2

例12电路图见图1-3。R=149,L=2H,C=0.32F。初始状态:电感电流为零,电容 电压为0.5V。t=0时刻接入1伏的电压。求0<t<15秒时i(t)和vt的值,并且画出电流与 电容电压的关系曲线。 根据基尔霍夫定律[RLC]有 di Ri+L-+V。=V 1.4 0.32F 令 得 以上方程式在M文件 electsys.m中定义如下: function xdot= Felectsys(t,x),%状态导数 V=1; %阶跃输入 R=14:L=2:C=0.32 dot=x(2)C, 1/L *(V-x(1)-R*X(2))1 下面的M-文件使用函数ode23对系统进行仿真 t0=0 tfinal=15 时间间隔0~15秒 x0=[0.5,0 %初值条件x1=0.5,x2=0 tol=0.001; %精度 %若非零值则打印出每一步的计算值 It, x=ode23('electsys', to, tfinal, xO, tol, trace); subplot(211),plot(t, x); of a rlc series circuit) xlabel(Time-sec'") text(8, 1.05, Capacitor voltage,), text(8, 0.05, current) vc=x(:,1),=x(,2) subplot(2 12),plot( vc, 1); title('current versus capacitor voltage") xlabel(Capacitor Voltage,) ylabel(current') subplot(lll) 仿真结果见图1-4 3

3 例 1.2 电路图见图 1-3。R=1.4Ω，L=2H，C=0.32F。初始状态：电感电流为零，电容 电压为 0.5V。ｔ＝０时刻接入１伏的电压。求 0<t<15 秒时 i(t)和 v(t)的值，并且画出电流与 电容电压的关系曲线。 根据基尔霍夫定律[RLC]有： v0 vs dt di Ri L   和 dt dv i C 0  令 x v , x i 1  0 2  得 1 2 x C 1 x  (v x Rx ) L 1 x2  s  1 2  以上方程式在 M-文件 electsys.m 中定义如下： function xdot=electsys(t, x); %状态导数 V=1; %阶跃输入 R=1.4; L=2; C=0.32; xdot=[x(2)/C;1/L*(V-x(1)-R*x(2))]; 下面的 M-文件使用函数 ode23 对系统进行仿真。 t0=0;tfinal=15; ﹡时间间隔 0～15 秒 x0=[0.5,0]; ％初值条件 x1=0.5,x2=0 tol=0.001; ％精度 trace=0; ％若非零值则打印出每一步的计算值 [t, x]=ode23('electsys',t0,tfinal,x0,tol,trace); subplot(211), plot(t, x); title('Time response of a RLC series circuit') xlabel('Time-sec') text(8,1.05,'Capacitor voltage'),text(8,0.05,'current') vc=x(:,1), i=x(:,2); subplot(212),plot(vc, i); title('current versus capacitor voltage') xlabel('Capacitor Voltage') ylabel('current') subplot(111) 仿真结果见图 1-4。 图 1-3

Time response of a RLC series circuit 15 0.5 Time-sec current versus capacitor voltage 0.15 0.1 05 -0.05 0.50.60.7 0.8 0.9 1.1 1.2 1.3 图1 113非线性系统 在变量的一定范围内大多数物理系统是线性的。然而,随着变量范围的延伸,所有的系 统最终都是非线性的。对于非线性系统,叠加原理不再成立。Ode23和ode45可以求解非线 性微分方程,见例1-3 例13图1-5的单摆,一长L米的无重量绳悬挂重为W=mg千克 的物体阻尼系数为B千克/米/秒。通常,该系统用线性微分方程近似描 述,实际上它是非线性的。 如果θ为绳的倾角,物体的速度将为Lθ,那么离心力为 F=-Wsin0-ble 根据牛顿定理有 图1 F=mle 联结上两式,得 ml 0+bl0+ wsin a=0 令x1=6,x2=6得

4 0 5 10 15 -0.5 0 0.5 1 1.5 Time response of a RLC series circuit Time-sec Capacitor voltage current 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 -0.1 -0.05 0 0.05 0.1 0.15 current versus capacitor voltage Capacitor Voltage current 图 1-4 1.1.3 非线性系统 在变量的一定范围内大多数物理系统是线性的。然而，随着变量范围的延伸，所有的系 统最终都是非线性的。对于非线性系统，叠加原理不再成立。Ode23 和 ode45 可以求解非线 性微分方程，见例 1-3。 例 1.3 图 1-5 的单摆，一长Ｌ米的无重量绳悬挂重为 W=mg 千克 的物体,阻尼系数为 B 千克/米/秒。通常，该系统用线性微分方程近似描 述，实际上它是非线性的。 如果θ为绳的倾角，物体的速度将为 Lθ，那么离心力为     W BL F sin - 根据牛顿定理有  F mL 联结上两式，得 mL  BL  Wsin   0   令      1 x2 x , 得

上面方程式在M-文件 pendulum.m中定义如下 function dot=pe %状态导数 W=2;L=6,B=0.02;g=9.81;m=W/g; dot=x(2),-B/m*x(2)-W/(m)*sin(x(1) 下面的M-文件对系统进行仿真 0=0: tfinal=s %时间间隔0~5秒 初始条件x1=1,x2=0 tol=0 0001: trace=0 %精度,不打印每一步的计算值 It, x]=ode23( pendulum, to, tfinal, xo, tol, trace); subplot(211),plot(t, x) title( Time response of a rigid pendulum,),xlabel(Time-sed text(3. 2, 3.5, Velocity,), text(3.2, 1.2, Angle-Rad') th=x(,1),w=x(;2) subplot(212), plot(th, w) lane plot of pendulum) xlabel('Position-Rad), ylabel(,Angular velocity) subplot(lll 仿真结果见图1-6 Time response of a rigid pendulum 0.5 Phase plane plot of pendulum 1-0.80.60.40.200.20.40.60.81 Position-Rad 图1-6

5 1 2 2 2 1 sinx mL W x m B x  x , x    上面方程式在 M-文件 pendulum.m 中定义如下： function xdot=pendulum(t, x); ％状态导数 W=2; L=.6; B=0.02; g=9.81; m=W/g; xdot=[x(2); -B/m*x(2)-W/(m*L)*sin(x(1))]; 下面的 M-文件对系统进行仿真。 t0=0; tfinal=5; ％时间间隔 0～5 秒 x0=[1,0]; ％初始条件 x1=1,x2=0 tol=0.0001; trace=0; ％精度，不打印每一步的计算值 [t, x]=ode23('pendulum',t0,tfinal,x0,tol,trace); subplot(211),plot(t, x) title('Time response of a rigid pendulum'),xlabel('Time-sec') text(3.2,3.5,'Velocity'),text(3.2,1.2,'Angle-Rad') th=x(:,1);w=x(:,2); subplot(212),plot(th, w) title('Phase plane plot of pendulum') xlabel('Position-Rad'),ylabel('Angular velocity') subplot(111) 仿真结果见图 1-6。 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 -4 -2 0 2 4 Time response of a rigid pendulum Time-sec Velocity Angle-Rad -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 -4 -2 0 2 4 Phase plane plot of pendulum Position-Rad Angular velocity 图 1-6

1.14线性化 在小信号条件下,非线性系统可以线性化。例1.3中单摆的运动方程在初始角较小的 情况下,可被线性化。设0=0+△0,单摆运动方程如下: mL(6+△6)+BL(0+△6)+wsin(。+△)=0 当△0(《∞时,有sin△0≈0,cosΔ0≈1,又因为0o较小,有sin0o≈Δ0,那么有线 性方程 mL△e+BLA+W△6=0 当Δθ很小时,线性方程(12)的解与(1.1)一致。 1.2传递函数 线性时不变系统的传递函数定义为零初始条件下输出量的拉普拉斯变换与输入量的拉 普拉斯变换像函数之比。尽管传递函数只能用于线性系统,但它比微分方程提供更为直观的 信息。令传递函数的分母多项式等于零,便得到特征方程。特征方程的根是系统的极点,分 子多项式的根是系统的零点。那么传递函数便可由常数项与系统的零、极点确定。常数项, 通常记作k,是系统的增益。利用传递函数,我们可以方便的研究系统参数的改变对响应的 影响。通过拉普拉斯反变换可得到系统在时域的响应。通常需要用有理函数的部分分式展开 在这部分举几个例子介绍 MATLLAB中求特征多项式的根,求传递函数的零、极点, 部分分式展开以及已知零、极点求传递函数等函数的功能 121多项式的根和特征多项式 如果P是包含多项式系数的行向量,roos(P)得到一个列向量,其元素为多项式的根。 如果r是包含多项式根的一个行例列向量,poly(r)得到一个行向量,其元素为多项式的系数 例14求多项式 s6+9s5+3125s4+6125s3+67.752+1475+15的根 多项式系数以降幂次序排列在一行向量中。用 roots求根。 >>P=[19312561.2567.7514.7515 >>r=roots(P) 多项式的根从列向量r中得到 -40000 -3.0000 1.0000+2.00001 -0.0000+500001 -0.0000-5.0000i 例1.5多项式的根为-1,-2,-3±14。求多项式方程

6 1.1.4 线性化 在小信号条件下，非线性系统可以线性化。例１.３中单摆的运动方程在初始角较小的 情况下，可被线性化。设θ=θ+Δθ，单摆运动方程如下： mL(0  )  BL(0  )  Wsin(0  )  0     (1.1) 当Δθ〈〈∞时，有 sinΔθ≈0,cosΔθ≈1,又因为θ0较小，有 sinθ0≈Δθ，那么有线 性方程 mL   BL  W   0   (1.2) 当Δθ很小时，线性方程（1.2）的解与（1.1）一致。 1.2 传递函数 线性时不变系统的传递函数定义为零初始条件下输出量的拉普拉斯变换与输入量的拉 普拉斯变换像函数之比。尽管传递函数只能用于线性系统，但它比微分方程提供更为直观的 信息。令传递函数的分母多项式等于零，便得到特征方程。特征方程的根是系统的极点，分 子多项式的根是系统的零点。那么传递函数便可由常数项与系统的零、极点确定。常数项， 通常记作 k，是系统的增益。利用传递函数，我们可以方便的研究系统参数的改变对响应的 影响。通过拉普拉斯反变换可得到系统在时域的响应。通常需要用有理函数的部分分式展开。 在这部分举几个例子介绍 MATLLAB 中求特征多项式的根，求传递函数的零、极点， 部分分式展开以及已知零、极点求传递函数等函数的功能。 1.2.1 多项式的根和特征多项式 如果 P 是包含多项式系数的行向量，roots(P)得到一个列向量，其元素为多项式的根。 如果 r 是包含多项式根的一个行/列向量，poly(r)得到一个行向量，其元素为多项式的系数。 例 1.4 求多项式 s 6+9s 5+31.25s 4+61.25s 3+67.75s 2+14.75s+15 的根。 多项式系数以降幂次序排列在一行向量中。用 roots 求根。 >>P=[1 9 31.25 61.25 67.75 14.75 15]; >>r=roots(P) 多项式的根从列向量 r 中得到 r = -4.0000 -3.0000 -1.0000 + 2.0000i -1.0000 - 2.0000i -0.0000 + 5.0000i -0.0000 - 5.0000i 例 1.5 多项式的根为-1,-2,-3±j4。求多项式方程

为了输入复数,必须首先建立虚数单位。然后在行列向量中输入根。使用poly得到多 项式方程 >>=sqrt(-1) >>r=[-1;2-3+4*-3-4* > p=poly(r) 多项式的系数从行向量中得到 因此,多项式方程为 s4+9s3+45s2+875+50=0 例16求下列矩阵的特征方程的根 A=-6-116 -6-115 用ploy求矩阵的特征方程的根。用 roots求方程的根 >>A=[01-1;-6-116;-6-115] >>P=poly (A) >>Froots(P) 结果如下 10000600001100006.0000 30000 20000 -1.0000 122传递函数的零点和极点 ·函数求传递函数的零点,极点和增益。 例17求下列传递函数的零点,极点和增益。 s3+1ls2+30s H(s)= s+9s+45s2+87s+50 >>num=[111300 >>den={19458750]; >>[z, p, k=tf2zp(num, den

7 为了输入复数，必须首先建立虚数单位。然后在行/列向量中输入根。使用 poly 得到多 项式方程。 >> i=sqrt(-1); >> r=[-1;-2;-3+4*i;-3-4*i]; >> p=poly(r) 多项式的系数从行向量中得到 p = 1 9 45 87 50 因此，多项式方程为 s 4+9s 3+45s 2+87s+50=0 例 1.6 求下列矩阵的特征方程的根        - 6 -11 5 - 6 -11 6 0 1 -1 A 用 ploy 求矩阵的特征方程的根。用 roots 求方程的根。 >> A=[0 1 -1;-6 -11 6;-6 -11 5]; >> P=poly(A) >> r=roots(P) 结果如下： P= 1.0000 6.0000 11.0000 6.0000 r = -3.0000 -2.0000 -1.0000 1.2.2 传递函数的零点和极点  函数 tf2zp 求传递函数的零点，极点和增益。 例 1.7 求下列传递函数的零点，极点和增益。 s 9s 45s 87s 50 s 11s 30s H(s) 4 3 2 3 2        >> num=[1 11 30 0]; >> den=[1 9 45 87 50]； >> [z, p, k]=tf2zp(num, den)

0 -3.0000-4.00001 因而有 s(S+5)(s+6) Hs)(s+1)s+2s+3+j4)(s+3-j4) ·函数zp2tf根据给定零点,极点和增益求传递函数。 例1.8系统的零点为-6,-5,0,极点为-3±14,-2,-1增益为1。求其传递函数。 -2;-1] >>num, den=zp2tf(z, p, k) 上面程序的结果为 因此,传递函数为 s3+1ls2+30s H(s) s4+9s3+45s2+87s+50 123部分分式展开 函数[r,p,k]F= residue(b,a,对两个多项式的比进行部分分式展开,如 Q(s)a (1.3) 向量b,a以s的降幂顺序排列多项式的系数。部分分式展开后余数送入列向量r,极点 送入列向量p,常数项送入k 例1.9对F(s)进行部分分式展开

8 z = 0 -6.0000 -5.0000 p = -3.0000 + 4.0000i -3.0000 - 4.0000i -2.0000 -1.0000 因而有 (s 1)(s 2)(s 3 j4)(s 3 j4) s(s 5)(s 6) H(s)           函数 zp2tf 根据给定零点，极点和增益求传递函数。 例 1.8 系统的零点为-6,-5,0,极点为-3±j4,-2,-1,增益为 1。求其传递函数。 >> z=[-6;-5;0];k=1; >> i=sqrt(-1); >> p=[-3+4*i;-3-4*i;-2;-1]; >> [num, den]=zp2tf(z, p, k) 上面程序的结果为 num = 0 1 11 30 0 den = 1 9 45 87 50 因此，传递函数为 s 9s 45s 87s 50 s 11s 30s H(s) 4 3 2 3 2        1.2.3 部分分式展开 函数[r, p, k]=residue(b, a)，对两个多项式的比进行部分分式展开，如 a s a s a s a b s b s b s b Q(s) P(s) 1 n 1 n 1 n n 1 0 m 1 m 1 m m               (1.3) 向量 b, a 以 s 的降幂顺序排列多项式的系数。部分分式展开后余数送入列向量 r，极点 送入列向量 p，常数项送入 k。 例 1.9 对 F(s)进行部分分式展开