现代控制理论实验报告

现代控制理论实验报告

二〇一六年五月

实验一 线性定常系统模型

一 实验目的

1. 掌握线性定常系统的状态空间表达式。学会在MATLAB中建立状态空间模型的方法。

2. 掌握传递函数与状态空间表达式之间相互转换的方法。学会用MATLAB实现不同模型之间的相互转换。

3. 熟悉系统的连接。学会用MATLAB确定整个系统的状态空间表达式和传递函数。

4. 掌握状态空间表达式的相似变换。掌握将状态空间表达式转换为对角标准型、约当标准型、能控标准型和能观测标准型的方法。学会用MATLAB进行线性变换。

二 实验内容

1. 已知系统的传递函数

G(s)?4 2s(s?1)(s?3)(1）建立系统的TF或ZPK模型。

（2）将给定传递函数用函数ss( )转换为状态空间表达式。再将得到的状态空间表达式用函数tf( )转换为传递函数，并与原传递函数进行比较。

（3）将给定传递函数用函数jordants( )转换为对角标准型或约当标准型。再将得到的对角标准型或约当标准型用函数tf( )转换为传递函数，并与原传递函数进行比较。

（4）将给定传递函数用函数ctrlts( )转换为能控标准型和能观测标准型。再将得到的能控标准型和能观测标准型用函数tf( )转换为传递函数，并与原传递函数进行比较。

2. 已知系统的传递函数

1??0?0??x???x??1?u y??11?x

?5?6????

（1）建立给定系统的状态空间模型。用函数eig( ) 求出系统特征值。用函数tf( ) 和zpk( )将这些状态空间表达式转换为传递函数，记录得到的传递函数和它的零极点。比较系统的特征值和极点是否一致，为什么?

（2）用函数canon( )将给定状态空间表达式转换为对角标准型。用函数eig( )求出系统特征值。比较这些特征值和（1）中的特征值是否一致，为什么? 再用函数tf( )和zpk( )将对角标准型或约当标准型转换为传递函数。比较这些传递函数和（1）中的传递函数是否一致，为什么?

（3）用函数ctrlss( )将给定的状态空间表达式转换为能控标准型和能观测标准型。用函数eig( )求系统的特征值。比较这些特征值和（1）中的特征值是否一致，为什么?再用函数tf( )将它们转换为传递函数。比较这些传递函数和（1）中的传递函数是否一致，为什么?

三 实验结果

1. 已知系统的传递函数 G(s)?

>> num=4; >> den=[1 5 7 3 0]; >> G=tf(num,den)

Transfer function: 4 ------------------------- s^4 + 5 s^3 + 7 s^2 + 3 s >> z=[]; >> p=[0 -1 -1 -3]; >> k=4; >> G=zpk(z,p,k)

Zero/pole/gain: 4 --------------- s (s+1)^2 (s+3)

4 2s(s?1)(s?3)

>> Gss=ss(G) a =

x1 x2 x3 x4 x1 0 1 0 0 x2 0 -1 1 0 x3 0 0 -1 1 x4 0 0 0 -3 b = u1 x1 0 x2 0 x3 0 x4 2 c =

x1 x2 x3 x4 y1 2 0 0 0 d = u1 y1 0

Gtf=tf(Gss)

Transfer function: 4 ------------------------- s^4 + 5 s^3 + 7 s^2 + 3 s

>> Gzpk=zpk(Gtf)

Zero/pole/gain:

4 --------------- s (s+3) (s+1)^2

1??0?0????2. xx?u y??11?x ?????5?6??1?

>> A=[0 1; -5 -6 ];B=[0;1]; >> C=[1 1];D=0; >> G=ss(A,B,C,D) a =

x1 x2 x1 0 1 x2 -5 -6 b = u1 x1 0 x2 1 c =

x1 x2 y1 1 1 d = u1 y1 0

Continuous-time model.

>> a=[0 1 ;-5 -6] a =

>> A=[0 0 -1;1 0 -3;0 1 -3]; B=[1;1;0]; C=[0 1 -2]; [Ac Bc Cc Tc Kc]=ctrbf(A,B,C) Ac =

-1.0000 0.0000 -0.0000 -2.1213 -2.5000 0.8660 -1.2247 -2.5981 0.5000 Bc =

0 0 -1.4142 Cc =

1.7321 1.2247 -0.7071 Tc =

-0.5774 0.5774 -0.5774 0.4082 -0.4082 -0.8165 -0.7071 -0.7071 0 Kc =

1 1 0 2.

unction [Ak, Bk, Ck ,Tk]=kalmdec(A,B,C) [Ac,Bc,Cc,Tc,Kc]=ctrbf(A,B,C); nc=rank(ctrb(A,B));

%按能控能观测性分解 n=length(A); ic=n-nc+1:n;

[Ao1,Bo1,Co1,To1,Ko1]=obsvf(Ac(ic,ic),Bc(ic),Cc(ic)); if nc

[Ao2,Bo2,Co2,To2,Ko2]=obsvf(Ac(inc,inc),Bc(inc),Cc(inc)); end

[m1,n1]=size(To1); [m2,n2]=size(To2); To=[To2,zeros(m2,n1);

zeros(m1,n2),To1]; T=To*Tc;

n1=rank(obsv(Ac(ic,ic),Cc(ic))); n2=rank(obsv(Ac(inc,inc),Cc(inc)));

K=[zeros(1,n-nc-n2),ones(1,n2),2*ones(1,nc-n1),3*ones(1,n1)]; Ak1=T*A*inv(T);Bk1=T*B;Ck1=C*inv(T); Ak=rot90(Ak1,2); Bk=rot90(Bk1,2); Ck=rot90(Ck1,2); Tk=rot90(T,2);

>> A=[0 0 -1;1 0 -3;0 1 -3]; B=[1;1;0]; C=[0 1 -2]; [Ak Bk Ck Tk]=kalmdec(A,B,C) Ac =

-1.0000 0 1.2247 3.4641 -1.0000 -2.1213 -0.0000 0.0000 -1.0000 Bc = -0.7071 -1.2247 0 Cc =

-1.4142 -0.0000 1.7321 Tc =

0.7071 -0.0000 -0.7071 -0.4082 -0.8165 -0.4082

-0.5774 0.5774 -0.5774

实验四 稳定性

一 实验目的

掌握系统稳定性的概念。学会使用MATLAB确定线性定常系统和非线性定常系统的稳定性。

二 实验内容

1. 已知线性系统

?10???? x?x 01??（1）用函数eig( )，pole( )和zpkdata( )求出系统的特征值和极点。用函数pzmap( )绘制系统的零点和极点。确定系统的稳定性。

（2）任意给定对称正定矩阵Q，用函数lyap( )求解Lyaponov方程，确定系统的稳定性。与（1）的结果进行比较。

?0?（3）令B???，C??00? ，D?0，任意给定初始状态。用函数initial( )求出

?0?系统的零输入响应，并绘制相应的状态响应曲线。说明稳定系统的状态响应曲线与不稳定系统的状态响应曲线的区别。

?1?（4）令B???，C??11? ，D?0，初始状态为零。 用函数step( )求出系统

??1?在单位阶跃信号作用下的状态响应和输出响应, 并绘制相应的曲线。分析系统的状态稳定和输出稳定是否一致。

2. 已知非线性系统

??1??3x1?x2?x3 ??x2?x1?x2?x2编制相应的程序，用克拉索夫斯基法确定系统在原点处的稳定性。

三 实验结果

1.

(1) A=[1 0 ;0 1]; >> eig(A) ans = 1 1

（2）A=[0 1;-2 -3];Q=eye(2);P=lyap(A,Q) P =

1.0000 -0.5000 -0.5000 0.5000

>> s=posdef(P) s =

matrix is positive definite matrix

可知系统是稳定的。

(2) function s=posdef(P) %判定矩阵的正定性 r=length(P); for i=1:r

pp(i)=det(P(1:i,1:i)); end

k=find(pp<0);j=find(pp==0); if isempty(j)&isempty(k)

109876x 107Impulse Responsewith PID control:k=1,kd=1,ki=1Amplitude54321000.511.522.5Time (sec)33.544.55

初始PID参数摆角状态图

系统响应是不稳定的，不能满足要求，需要调整参数KP，KD和KI，直到获得满意的控制结果。首先增加比例系数KP，观察它对响应的影响，取KP=100，kd=1.系统响应如 下： num =

2.3566 0 0 den =

1.0000 0.0883 -27.8285 -2.3094 0 kd = 1 k = 100 ki = 1 numc =

2.3566 0 0 0 denc =

1.0000 2.4449 207.8268 0.0471 0 0

Impulse Responsewith PID control k=100,ki=1,kd=10.150.10.05Amplitude0-0.05-0.1-0.15-0.200.511.522.5Time (sec)33.544.55调节PID参数后摆角输出图

系统稳定时间约为4秒，满足要求。由于此时稳态误差为0，所以不需要改变积分环节（你可以改变积分系数，观察系统响应如何变化）；系统响应的超调量比较大，为了减小超调，增加微分系数KD，取KD=20，响应结果和响应曲线如下：

num =

2.3566 0 0 den =

1.0000 0.0883 -27.8285 -2.3094 0 kd = 20 k = 100 ki = 1 numc =

2.3566 0 0 0 denc =

1.0000 47.2194 207.8268 0.0471 0 0

Impulse Responsewith PID control k=100,ki=20,kd=10.10.090.080.070.06Amplitude0.050.040.030.020.01000.511.522.5Time (sec)33.544.55

微调PID参数后摆角状态曲线

s='matrix is positive definite matrix'; elseif isempty(k)

s='matrix is half positive definite matrix'; else

s='non-definite matrix'; end

>> A=[0 1;-2 -3];Q=eye(2);P=lyap(A,Q) s=posdef(P) 返回 P =

1.0000 -0.5000 -0.5000 0.5000 s =

matrix is positive definite matrix 2. %ex43 syms x1 x2; x=[x1 x2]; f1=-3*x1+x2; f2=x1-x2-x2^3;

F1=[jacobian(f1,x1) jacobian(f1,x2);jacobian(f2,x1) jacobian(f2,x2)]; F2=[jacobian(f1,x1) jacobian(f2,x1);jacobian(f1,x2) jacobian(f2,x2)] ; F=Ft+F2; n=length(x); for i=1:n

ff(i)=simple(det(F(1:i,1:i))); end ff

v=f1*f1+f2*f2

在命令窗中运行该程序得 ff =

[-6, 8+36*x2^2] v =

(-3*x1+x2)^2+(x1-x2-x2^3)^2 可见F负定, 且当x??时，

v(x)??，所以系统大范围渐近稳定。

实验五 极点配置和状态观测器

一 实验目的

1. 掌握状态反馈和输出反馈的概念及性质。

2. 掌握利用状态反馈进行极点配置的方法。学会用MATLAB求解状态反馈矩阵。

3. 掌握状态观测器的设计方法。学会用MATLAB设计状态观测器。 4. 熟悉分离定理，学会设计带有状态观测器的状态反馈系统。

二 实验内容

1. 已知系统

??300??1????020?x??1?u x???????00?1???1?? y??0.40.26670.333?x3

（1）求解系统的零点、极点和传递函数，并判断系统的能控性和能观测性。 （2）分别选取K=[0 3 0]，K=[1 3 2]，K=[0 16 /3 –1/3]为状态反馈矩阵，求解闭环系统的零点、极点和传递函数，判断闭环系统的能控性和能观测性。它们是否发生改变?为什么？

（3）任选三个输出反馈矩阵，求解闭环系统的零点、极点和传递函数，并判断系统的能控性和能观测性。它们是否发生改变? 为什么？ 2. 已知系统

0??01?0??x??0?u ???001 x???????0?2?3???1?? y??100?x

（1）求解系统的极点。绘制系统的单位阶跃响应曲线，并确定系统的超调量和上升时间。

（2）求解状态反馈矩阵K，使闭环系统的极点为?3和?13?j。求解状态反22馈系统的传递函数。绘制该闭环系统的单位阶跃响应曲线，并确定系统的超调量和上升时间。与原系统比较, 性能是否改善？

（3）设计一个全维观测器，使观测器的极点为-5，-5，-5。仿真状态观测器观测到的状态。

（4）设计一个降维观测器，使观测器的极点为-5。

（5）建立带全维状态观测器的状态反馈系统的状态空间表达式。求解带全维状态观测器的状态反馈系统的极点，是否是状态反馈系统和观测器的极点的组合？为什么？求解该闭环系统的传递函数，与状态反馈系统的传递函数是否一致？为什么？绘制该闭环系统的单位阶跃响应曲线，并确定系统的超调量和上升时间。与状态反馈系统的单位阶跃响应曲线比较，验证两种反馈是否等价。

三 实验结果

1. >> A=[-1 0 0;0 0 1;0 -3 1];B=[0 0 1]';Uc=ctrb(A,B);m=rank(Uc) m =

可见系统是不完全能控的，而且其中有两个极点可以任意配置。在命令窗中调用此函数

>> K=bfjdpz(A,B,C) 返回 Af =

-1 0 0 0 0 -1 0 3 1 Bf = 0 0 -1 Cf =

-1 1 -1 Tf =

-1 0 0 0 1 0 0 0 -1 输入能配置的极点组成的向量

由上述能控性分解可见，极点-1是不能控的，另外两个极点是能控的。可判

断利用状态反馈第一组极点可以配置，而第二组极点不能配置。在命令窗中输入第一组极点中能控的两个极点，即 [-2 -2] 返回 Pc =

-2 -2 K =

-1 1 5

2.

>> A=[-1 1;0 -2];B=[0 1]';C=[2 0];p=[-10 -10];L=(acker(A',C',p))’ 返回 L = 8.5000 32.0000

利用函数simobsv( )可以仿真状态观测器所观测到的状态。该函数的调用格式为

[Xo,X,t]?simobsv(G,L)

其中G为系统的状态空间模型，L为观测器的反馈矩阵。Xo为重构状态的阶跃响应，x为原系统状态的阶跃响应，t为函数自动选择的时间向量。该函数的程序如下 function [xo,x,t]=simobsv(G,L) %仿真观测器的状态 [y,t,x]=step(G);

G=ss(G);A=G.a;B=G.b;C=G.c;D=G.d; [y1,xo1]=step((A-L*C),B,C,D,1,t); [y2,xo2]=lsim((A-L*C),L,C,D,y,t); xo=xo1+xo2;

现在仿真上面的求得的状态观测器的状态。运行下列命令 >> figure('pos',[50 50 200 150],'color','w')； [xo,x,t]=simobsv(ss(A,B,C,0),L);plot(t,x,'-k',t,xo,':r') 返回原系统和状态观测器的单位阶跃响应曲线，如图1。

图1 单位阶跃响应曲线 图2 反馈矩阵响应曲线 再选观测器的反馈矩阵为L??面的效果好。

??2?，运行上述命令得到图2。可见此时的估计不如上??20?实验六 倒立摆

一 实验目的

倒立摆是研究控制理论的一种典型实验装置，具有成本低，结构简单，物理参数和结构易于调整的优点，是一个具有高阶次、不稳定、多变量、非线性和强耦合特性的不稳定系统。在控制过程中，他能有效反应鲁棒性，随动性，等控制中的关键问题。研究倒立摆可以让我们掌握和理解控制科学中很多的理论问题。

二 实验内容

倒立摆系统的构成

图1 倒立摆系统的组成框图

如图1所示为倒立摆的结构图。系统包括计算机、运动控制卡、伺服机构、倒立摆本体和光电码盘几大部分，组成了一个闭环系统。光电码盘1将小车的位移、速度信号反馈给伺服驱动器和运动控制卡，摆杆的位置、速度信号由光电码盘2反馈回控制卡。计算机从运动控制卡中读取实时数据，确定控制决策（小车向哪个方向移动、移动速度、加速度等），并由运动控制卡来实现该控制决策，产生相应的控制量，使电机转动，带动小车运动，保持摆杆平衡。

单级倒立摆数学模型的建立

单级倒立摆模型：

图1 单级倒立摆模型示意图

那我们在本实验中定义如下变量：

M 小车质量 （本实验系统 0.5 Kg） m 摆杆质量 （本实验系统 0.2 Kg） b 小车摩擦系数 （本实验系统 0.1 N/m/sec） l 摆杆转动轴心到杆质心的长度（0.3 m） I 摆杆惯量 （0.006 kg*m*m） F 加在小车上的力 x 小车位置

φ 摆杆与垂直向上方向的夹角

θ 摆杆与垂直向下方向的夹角（考虑到摆杆初始位置为竖直向下）

应用Newton方法来建立系统的动力学方程过程如下： 分析小车水平方向所受的合力，可以得到以下方程：

??F?bx??N M?x由摆杆水平方向的受力进行分析可以得到下面等式：

N?md2dt2(x?lsin?)

??cos??ml??2sin? ??ml?x即 N?m?把这个等式代入上式中，就得到系统的第一个运动方程：

??cos??ml??2sin??F （1） ??bx??ml?(M?m)?x为了推出系统的第二个运动方程，我们对摆杆垂直方向上的合力进行分析，可以得到下面方程：

d2P?mg?m2(lcos?)dt??sin??ml??2cos?即：P?mg??ml?力矩平衡方程如下：

?? ?Plsin??Nlcos??I?

注意：此方程中力矩的方向，由于?????,cos???cos?,sin???sin?，故等式前面有负号。

2合并这两个方程，约去P和N，由I?1得到第二个运动方程： 3ml43???mglsin???ml??cos? （2） ml2?x设?????（?是摆杆与垂直向上方向之间的夹角），假设?与1（单位是弧

sin????，(度）相比很小，即?《1，则可以进行近似处理：cos???1，

d?2)?0。dt用u来代表被控对象的输入力F，线性化后两个运动方程如下：

4???g?????x?3l? ??????bx??ml??ux?(M?m)?传递函数

对方程组(3)进行拉普拉斯变换，得到

?422?l?(s)s?g?(s)?X(s)s ?3?(M?m)X(s)s2?bX(s)s?ml?(s)s2?U(s)?注意：推导传递函数时假设初始条件为0。

由于输出为角度?，求解方程组（4）的第一个方程，可以得到

4gX(s)?[l?2]?(s)

3s把上式代入方程组（4）的第二个方程，得到

?(I?ml2)(M?m)??ml??(I?ml2)g?g?2?2??(s)s?ml?(s)s2?U(s) ??(s)s?b?s?s??ml整理后得到传递函数：

ml2sqs4?43?(s)?U(s)bmlq2

s3?(M?m)mgl2bmgls?sqq其中 q?[(M?m)(I?ml2)?(ml)2] （1）观察系统Matlab 仿真和开环响应

（2）观察稳定性与可控性分析

（3）PID算法：调整参数看仿真结果

三 实验结果

在Matlab中，拉普拉斯变换后得到的传递函数可以通过计算并输入分子和分母矩阵来实现。求系统传递函数的m-文件内容如下： M = 1.096; m = 0.109; b = 0.1; I= 0.0034; g = 9.8; l = 0.25;

q = (M+m)*(I+m*l^2)-(m*l)^2; %simplifies input

num = [m*l/q 0]

den = [1 b*(I+m*l^2)/q -(M+m)*m*g*l/q -b*m*g*l/q] t = 0 : 0.01 : 5;

impulse ( num , den , t )

本文来源：https://www.bwwdw.com/article/ywm2.html