控制系统计算机辅助设计 CAD 实验

控制系统计算机辅助

设计实验

一、实验教学目标与基本要求

上机实验是本课程重要的实践教学环节。实验的目的不仅仅是验证理论知

识，更重要的是通过上机加强学生的实验手段与实践技能，掌握应用MATLAB/Simulink 求解控制问题的方法，培养学生分析问题、解决问题、应用知识的能力和创新精神，全面提高学生的综合素质。

通过对MATLAB/Simulink进行求解，基本掌握常见控制问题的求解方法与命令调用，更深入地认识和了解MATLAB语言的强大的计算功能与其在控制领域的应用优势。

二、题目及解答

第一部分：MATLAB 必备基础知识、控制系统模型与转换、线性

控制系统的计算机辅助分析：

编写rossler.m文件：

function dx=rossler(t,x)

dx=[-x(2)-x(3);x(1)+0.2*x(2);0.2+(x(1)-5.7)*x(3)];

主函数：

>> x0=[0;0;0];

[t,y]=ode45('rossler',[0,100],x0);

plot3(y(:,1),y(:,2),y(:,3)),grid

>> plot3(y(:,1),y(:,2),y(:,3)),grid

>> plot(y(:,1),y(:,2)),grid

所的图像：三维相轨迹（下图）

xoy平面投影

2.

编写c2exmobj.m文件：

function y=c2exmobj(x)

y=x(1)^2-2*x(1)+x(2);

编写c2exmcon.m文件：

function [c,ce]=c2exmcon(x)

ce=[];

c=[4*x(1)^2+x(2)^2-4];

主函数为：

A=[];B=[];Aeq=[];Beq=[];xm=[0;0];xM=[];x0=[0;0];

ff=optimset;ff.Tolx=1e-10;ff.TolFun=1e-20;

x=fmincon('c2exmobj',x0,A,B,Aeq,Beq,xm,xM,'c2exmcon',ff)

运行结果：

Warning: Trust-region-reflective method does not currently solve this type of problem,

using active-set (line search) instead.

> In fmincon at 422

Optimization terminated: first-order optimality measure less than options.TolFun and maximum constraint violation is less than options.TolCon.

Active inequalities (to within options.TolCon = 1e-006):

lower upper ineqlin ineqnonlin

2 1

x =

1.0000

3.

(a).程序代码为：

>> s=tf('s');

G=(s^3+4*s+2)/(s^3*(s^2+2)*((s^2+1)^3+2*s+5))

运行结果为：

Transfer function:

s^3 + 4 s + 2

------------------------------------------------------ s^11 + 5 s^9 + 9 s^7 + 2 s^6 + 12 s^5 + 4 s^4 + 12 s^3 (b). 程序代码为：

>> z=tf('z',0.1);

H=(z^2+0.568)/((z-1)*(z^2-0.2*z+0.99))

运行结果为：

Transfer function:

z^2 + 0.568

-----------------------------

z^3 - 1.2 z^2 + 1.19 z - 0.99

Sampling time: 0.1

4.

将方程两边进行拉式变换后可得传递函数，代码如下：

>> tf('s');

G=2/(s^3+13*s^2+4*s+5)

Transfer function:

2

----------------------

s^3 + 13 s^2 + 4 s + 5

转换为状态空间方程为：

>> G1=ss(G)

a =

x1 x2 x3

x1 -13 -1 -1.25

x2 4 0 0

x3 0 1 0

b =

u1

x1 0.5

x2 0

x3 0

c =

x1 x2 x3

y1 0 0 1

d =

u1

y1 0

Continuous-time model.

转换为零极点模型为：

>> G2=zpk(G)

Zero/pole/gain:

2

----------------------------------

(s+12.72) (s^2 + 0.2836s + 0.3932)

5.

如题可设采样时间为0.1s，代码如下：

>> z=tf('z',0.1);

H=(z+2)/(z^2+z+0.16)

运行结果为：

Transfer function:

z + 2

--------------

z^2 + z + 0.16

Sampling time: 0.1

6.

编写feedback.m函数：

function H=feedback(G1,G2,key)

if nargin==2; key=-1; end, H=G1/(sym(1)-key*G1*G2); H=simple(H); 代码如下：

syms J Kp Ki s;

G=(s+1)/(J*s^2+2*s+5);Gc=(Kp*s+Ki)/s;

GG=feedback(G*Gc,1);G1=collect(GG,s)

运行结果为：

G1 =(s+1)*(Kp*s+Ki)/(J*s^3+(Kp+2)*s^2+(Ki+Kp+5)*s+Ki)

7.

(a):代码如下：

>> s=tf('s');G=(211.87*s+317.64)/(s+20)/(s+94.34)/(s+0.1684); Gc=(169.6*s+400)/s/(s+4);H=1/(0.01*s+1);GG=feedback(G*Gc,H)

运行结果为：

Transfer function:

359.3 s^3 + 3.732e004 s^2 + 1.399e005 s + 127056

----------------------------------------------------------------------------------

0.01 s^6 + 2.185 s^5 + 142.1 s^4 + 2444 s^3 + 4.389e004 s^2 + 1.399e005 s + 127056 状态方程如下：

>> ss(GG)

a =

x1 x2 x3 x4 x5 x6

x1 -218.5 -111.1 -29.83 -16.74 -6.671 -3.029

x2 128 0 0 0 0 0

x3 0 64 0 0 0 0

x4 0 0 32 0 0 0

x5 0 0 0 8 0 0

x6 0 0 0 0 2 0

b =

u1

x1 4

x2 0

x3 0

x4 0

x5 0

x6 0

c =

x1 x2 x3 x4 x5 x6

y1 0 0 1.097 3.559 1.668 0.7573

d =

u1

y1 0

Continuous-time model.

零极点模型转换如下：

>> zpk(GG)

Zero/pole/gain:

35933.152 (s+100) (s+2.358) (s+1.499)

-----------------------------------------------------------------------

(s^2 + 3.667s + 3.501) (s^2 + 11.73s + 339.1) (s^2 + 203.1s + 1.07e004) (b).程序代码如下：

可先设系统采样周期为0.1s

>> z=tf('z',0.1);G=(35786.7*z^2+108444*z^3)/(4*z+1)/(20*z+1)/(1+74.04*z);

Gc=z/(1-z);H=z/(0.5-z);GG=feedback(G*Gc,H)

运行结果为：

Transfer function:

-108444 z^5 + 1.844e004 z^4 + 1.789e004 z^3

--------------------------------------------------------------------- 1.144e005 z^5 + 2.876e004 z^4 + 274.2 z^3 + 782.4 z^2 + 47.52 z + 0.5 Sampling time: 0.1

状态方程转换如下：

>> ss(GG)

a =

x1 x2 x3 x4 x5

x1 -0.2515 -0.00959 -0.1095 -0.05318 -0.01791

x2 0.25 0 0 0 0

x3 0 0.25 0 0 0

x4 0 0 0.125 0 0

x5 0 0 0 0.03125 0

b =

u1

x1 1

x2 0

x3 0

x4 0

x5 0

c =

x1 x2 x3 x4 x5

y1 0.3996 0.6349 0.1038 0.05043 0.01698

d =

u1

y1 -0.9482

Sampling time: 0.1

Discrete-time model.

零极点模型转化如下：

>> zpk(GG)

Zero/pole/gain:

-0.94821 z^3 (z-0.5) (z+0.33)

----------------------------------------------------------

(z+0.3035) (z+0.04438) (z+0.01355) (z^2 - 0.11z + 0.02396)

Sampling time: 0.1

8.

程序代码如下：

s=tf('s');c1=feedback(1/(s+1)*s/(s^2+2),(4*s+2)/(s+1)^2);c2=feedback(1/s^2,5 0);

c3=feedback(c1*c2,(s^2+2)/(s^3+14));G=3*c3;

G

运行结果为：

Transfer function:

3 s^6 + 6 s^5 + 3 s^

4 + 42 s^3 + 84 s^2 + 42 s

-----------------------------------------------------------------------------------------------------

s^10 + 3 s^9 + 55 s^8 + 175 s^7 + 300 s^6 + 1323 s^5 + 2656 s^4 + 3715 s^3 + 7732 s^2 + 5602 s + 1400

9.

程序代码如下：

>> s=tf('s');

G=(s+1)^2*(s^2+2*s+400)/(s+5)^2/(s^2+3*s+100)/(s^2+3*s+2500);

G1=c2d(G,0.01),

G2=c2d(G,0.1),

G3=c2d(G,1)

运行结果如下：

Transfer function:

4.716e-005 z^5 - 0.0001396 z^4 + 9.596e-005 z^3 + 8.18e-005 z^2 - 0.0001289 z + 4.355e-005

------------------------------------------------------------------------------------------

z^6 - 5.592 z^5 + 13.26 z^4 - 17.06 z^3 + 12.58 z^2 - 5.032 z + 0.8521 Sampling time: 0.01

Transfer function:

0.0003982 z^5 - 0.0003919 z^4 - 0.000336 z^3 + 0.0007842 z^2 - 0.000766 z + 0.0003214

----------------------------------------------------------------------------

---------

z^6 - 2.644 z^5 + 4.044 z^4 - 3.94 z^3 + 2.549 z^2 - 1.056 z + 0.2019 Sampling time: 0.1

Transfer function:

8.625e-005 z^5 - 4.48e-005 z^4 + 6.545e-006 z^3 + 1.211e-005 z^2 - 3.299e-006 z + 1.011e-007

--------------------------------------------------------------------------------------------

z^6 - 0.0419 z^5 - 0.07092 z^4 - 0.0004549 z^3 + 0.002495 z^2 - 3.347e-005 z + 1.125e-007

Sampling time: 1

绘制曲线如下：

>> subplot(221),step(G)

subplot(222),step(G1)

subplot(223),step(G2)

subplot(224),step(G3)

图像为：

10.

(a)程序代码如下：

>> s=tf('s');G=1/(s^3+2*s^2+s+2);eig(G)

运行结果如下：

ans =

-2.0000

-0.0000 + 1.0000i

-0.0000 - 1.0000i

可得系统为临界稳定

(b)程序代码如下：

>> num=1;den=[6 3 2 1 1];G=tf(num,den);eig(G)

运行结果为：

ans =

-0.4949 + 0.4356i

-0.4949 - 0.4356i

0.2449 + 0.5688i

0.2449 - 0.5688i

可得有一对共轭复根在右半平面，所以系统不稳定。

(c)程序代码如下：

>> num=1;den=[1 1 -3 -1 2];G=tf(num,den);eig(G)

运行结果为：

ans =

-2.0000

-1.0000

1.0000

1.0000

可得有两根在右半平面，故系统不稳定。

11.

均假设系统采样时间为0.1s

(a).程序代码如下：

>> num=[-3 2];den=[1 -0.2 -0.25 0.05];H=tf(num,den,'Ts',0.1);abs(eig(H)')

运行结果如下：

ans =

0.5000 0.5000 0.2000

可得系统稳定。

(b). 程序代码如下：

>> num=[3 -0.39 -0.09];den=[1 -1.7 1.04 0.268

0.024];H=tf(num,den,'Ts',0.1);abs(eig(H)')

运行结果如下：

ans =

1.1939 1.1939 0.1298 0.1298

可得系统不稳定。

12.

程序代码如下：

>> A=[-0.2 0.5 0 0 0;0 -0.5 1.6 0 0;0 0 -14.3 85.8 0;0 0 0 -33.3 100;0 0 0 0 -10]; B=[0;0;0;0;30];C=[];D=[];G=ss(A,B,C,D);eig(G)

运行结果如下：

ans =

-0.2000

-0.5000

-14.3000

-33.3000

-10.0000

可得系统稳定。

13.

程序代码如下：

>> A=[-5 2 0 0;0 -4 0 0;-3 2 -4 -1;-3 2 0 -4];X0=[1 2 0 1]';

syms t; Y=expm(A*t)*X0

运算结果为：（此为解析解）

Y=

-3*exp(-5*t)+4*exp(-4*t)

2*exp(-4*t)

18*exp(-4*t)-18*exp(-5*t)-18*t*exp(-4*t)+4*t^2*exp(-4*t)

10*exp(-4*t)-9*exp(-5*t)-8*t*exp(-4*t)

数值解代码如下：

f.m函数如下：

function dx=f(t,x);

A=[-5 2 0 0;0 -4 0 0;-3 2 -4 -1;-3 2 0 4];

dx=A*x

主函数如下：

>> x0=[1;2;0;1];[t,y]=ode45('f',[0,5],x0)

运行结果为：

dx =

1.0e+009 *

-0.0000

-0.0000

-0.2426

1.9411

14.

(a).程序代码如下：

>> s=tf('s');G=(s+6)*(s-6)/s/(s+3)/(s+4+4j)/(s+4-4j);

rlocus(G)

可得根轨迹图像为：

不存在K使得系统稳定。

(b).程序代码如下：

>> s=tf('s');G=(s^2+2*s+2)/(s^4+s^3+14*s^2+8*s);rlocus(G)

图像如下：

可得0

15.

代码如下：

采用1/3阶Padé近似拟合延迟项

代码如下：

编写pade_app.m文件：

function Gr=pade_app(c,r,k)

w=-c(r+2:r+k+1)';vv=[c(r+1:-1:1)';zeros(k-1-r,1)];

W=rot90(hankel(c(r+k:-1:r+1),vv));V=rot90(hankel(c(r:-1:1))); x=[1 (W\w)'];dred=x(k+1:-1:1)/x(k+1);

y=[c(1) x(2:r+1)*V'+c(2:r+1)];nred=y(r+1:-1:1)/x(k+1);

Gr=tf(nred,dred);

编写paderm.m文件：

function [n,d]=paderm(tau,r,k)

c(1)=1;for i=2:r+k+1,c(i)=-c(i-1)*tau/(i-1);end

Gr=pade_app(c,r,k);n=Gr.num{1}(k-r+1:end);d=Gr.den{1};

>> [n,d]=paderm(2,1,3);s=tf('s');G=(s-1)/(s+1)^5;rlocus(G*tf(n,d)) 运行结果：

由图得K∈(0,3.91)能够使得闭环系统稳定

16.

程序代码如下：

>>s=tf('s');G=8*(s+1)/(s^2*(s+15)*(s^2+6*s+10));

figure;bode(G),figure;nyquist(G),figure;nichols(G),grid

3个图如下：

由奈氏图不包围（-1，j0）点及奈氏稳定判据可以得出，该系统稳定。[gm,pm,wcg,wcp]=margin(G)

gm =

30.4686

pm =

4.2340

wcg =

1.5811

wcp =

0.2336

阶跃响应为：

可得系统稳定。

代码如下：假设采样周期为0.1s

由图可知，系统不稳定。

[gm,pm,wcg,wcp]=margin(G)

gm =

30.4686

pm =

4.2340

wcg =

1.5811

wcp =

0.2336

阶跃响应如下：

可见系统最后发散，故系统不稳定。

17.

程序代码如下：

>>s=tf('s');

G=100*(1+s/2.5)/(s*(1+s/0.5)*(1+s/50));

Gc=1000*(s+1)*(s+2.5)/((s+0.5)*(s+50));

nyquist(G*Gc),grid

由奈氏图可得，曲线不包围(-1,j0)点，而开环系统不含有不稳定极点，所以根据奈氏稳定判据闭环系统是稳定的。

用阶跃响应来验证得：

可得系统是稳定的。

第二部分：Simulink 在系统仿真中的应用、控制系统计算机辅助设计、控制工程中的仿真技术应用

Simulink仿真图如下：

由plot(tout,yout)

得到响应仿真图为：

微分方程数值解

代码如下：

>>x0=[1;0.5;0.5;0.2];

[t,y]=ode45('diff1',[0,10],x0);

plot(t,y(:,1))

画的图像为:

可得，与simulink搭成模块仿真的一样。

3.

仿真图搭建如下：

仿真结果为：

输出曲线及误差曲线

>> plot(tout,yout)

4．

由图可得以下方程：

5.

模块搭建如下：

提取系统模型为：

>> [A,B,C,D]=linmod('dianji');G=ss(A,B,C,D)

运行结果为：

a =

x1 x2 x3 x4 x5 x6 x7 x8 x9

x1 0 0 0 0 0 0 0 0 1.4 x2 0 -100 0 200 -0.88 0 11.76 0 0 x3 0 0 -100 0 0 0 0 0 1.4 x4 0 0 0 -100 0 0 0 0 0 x5 130 0 0 0 -100 0 0 0 0 x6 0 294.1 -29.41 0 0 -149.3 0 19.61 0 x7 0 0 0 100 -0.44 0 0 0 0 x8 0 100 -10 0 0 0 0 0 0 x9 -27.56 0 0 0 0 1.045e+004 0 0 -6.66 b =

u1

x1 0

x2 0

x3 0

x4 1

x5 0

x6 0

x7 0

x8 0

x9 0

c =

x1 x2 x3 x4 x5 x6 x7 x8 x9

y1 130 0 0 0 0 0 0 0 0

d =

u1

y1 0

Continuous-time model.

Continuous-time model.

阶跃响应和频率响应曲线为：

subplot(221),step(G),grid

subplot(222),bode(G),grid

subplot(223),nyquist(G),grid

subplot(224),nichols(G),grid

6.

>> s=tf('s');

G=210*(s+1.5)/((s+1.75)*(s+16)*(s+1.5+3*j)*(s+1.5-3*j));

Gc=52.5*(s+1.5)/(s+14.86);sisotool(G)

此为系统校正前的图像：

由：sisotool(G，Gc)得校正后的图像：

有以下代码可得：

>> [gm,pm,wg,wp]=margin(G)

gm =

4.8921

pm =

60.0634

wg =

7.9490

wp =

3.9199

稳定裕量为60.0634°，可得系统是稳定的。

而加入控制器后：

>> [gm,pm,wg,wp]=margin(G*Gc)

gm =

0.8090

pm =

-6.0615

wg =

17.1659

wp =

18.9029

可得稳定裕量为-6.0615°，系统俨然不稳定了。

7.

可得状态反馈矩阵和Riccati方程的解为

>>A=[0,1,0,0;0,0,1,0;-3,1,2,3;2,1,0,0];B=[1,0;2,1;3,2;4,3];Q=diag([1,2,3,4])

;R=eye(2);[K,S]=lqr(A,B,Q,R)

K =

-0.0978 1.2118 1.8767 0.7871

-3.8819 -0.4668 2.6713 1.0320

S =

5.4400 0.6152 -2.3163 0.0452

0.6152 1.8354 -0.0138 -0.7582

-2.3163 -0.0138 1.9214 -0.3859

0.0452 -0.7582 -0.3859 0.8540

得闭环系统的极点为：

>> eig(A-B*K)

ans =

-12.2563

-1.6786 + 0.9981i

-1.6786 - 0.9981i

-1.4627

闭环系统各个状态的曲线

>> step(ss(A-B*K,B,eye(4),zeros(4,2)))

8.

求原系统的极点和零点

>>

A=[-0.2,0.5,0,0,0;0,-0.5,1.6,0,0;0,0,-14.3,85.8,0;0,0,0,-33.3,100;0,0,0,0,-1 0];B=[0;0;0;0;30];C=[1,0,0,0,0];D=0;G=ss(A,B,C,D);pole(G),zero(G)

ans =

-0.2000

-0.5000

-14.3000

-33.3000

-10.0000

ans =

Empty matrix: 0-by-1

配置系统极点

>> p=[-1,-2,-3,-4,-5];K=place(A,B,p),eig(A-B*K)

K =

0.0004 0.0004 -0.0035 0.3946 -1.4433

ans =

-5.0000

-4.0000

-3.0000

-2.0000

-1.0000

9.

建立模型xiti9.mdl

运行optimPID，加载xiti9.mdl后点击Create File创建函数optpidfun_0.m function y=optpidfun_0(x)

opt=simset('OutputVariables','y');

assignin('base','Kp',x(1));

assignin('base','Ki',x(2));

assignin('base','Kd',x(3));

try

[t_time,x_state,y_out]=sim('pidctrl_model',[0,30.000000],opt);

catch, y_out=10000; end

y=y_out(end,1);

设置相应参数后，点击Optimize后得到如下结果

ctrl_pars =

0.5326

0.2294

0.3561

yy =

1.2446

即比例，积分，微分相应系数为：Kp =0.5326，Ki = 0.2294，Kd = 0.3561再代入模型xiti9.mdl进行仿真运行，得到曲线如下图

>>plot(tout,yout)有：