华电线性系统理论大作业

分数: ___________

任课教师签字：___________ 华北电力大学研究生结课作业

学年学期：2014-2015学年第一学期

课程名称：线性系统理论

学生姓名：

学号：

提交时间：2014年11月27日

目录

1.绪论 (1)

2.球杆系统分析与建模 (1)

2.1球杆模型简介 (1)

2.2拉格朗日法建模 (1)

2.3拉格朗日模型线性化及状态空间表达式求取 (4)

3. 系统稳定性分析 (5)

3.1有初始状态下求取系统响应曲线 (6)

3.3稳定性判断并求取零极点分布图 (7)

4.系统能控性判别 (8)

4.1代数判据 (8)

4.2模态判据 (8)

4.3可控性与可稳定性 (10)

5.系统极点配置 (10)

5.1极点配置方法 (10)

5.1.1状态反馈原理 (11)

5.1.2输出反馈原理 (11)

5.1.3PID配置极点原理 (12)

5.1.4三种反馈对比 (12)

5.2.用状态反馈进行极点配置 (12)

6.可观性分析及带状态反馈的状态观测器的设计 (16)

6.1能观性分析 (16)

6.1.1代数判据 (16)

6.1.2模态判据 (16)

6.3全维观测器原理 (17)

6.4全维状态观测器结构 (17)

6.5全维状态观测器设计 (18)

6.6全维状态观测器Simulink仿真 (18)

6.7全维状态观测器在干扰下的性能研究 (20)

7.总结 (22)

1.绪论

球杆系统是控制理论中很经典的一个模型，通常用来检验控制策略的效果，并且很多实际系统都可以近似抽象为球杆模型，因此，对球杆系统的研究很有意义，本文从球杆模型的拉格朗日法建模入手，对球杆系统稳定性，能控能观性等控制特性进行分析。

2.球杆系统分析与建模

2.1球杆模型简介

球杆系统由底座，直流伺服电机，光滑导轨，小球等组成，导轨在伺服电机的带动下转动，小球在自身重力的作用下沿着光滑的金属导轨自由滚动，球杆系统简图如下，其中x 是小球在导轨上相对于导轨中心的位移量，以导轨左侧为正，

α是导轨相对于水平线的倾斜角。

图1.球杆系统简图

2.2拉格朗日法建模

为了对球杆系统进行研究，我们先对其进行建模，一般来说，这种球杆系统，运用拉格朗日方程建立其数学模型比较方便，拉格朗日方程如下：

()t U q R

q V q T q T dt d +??-??-=??-???? ????'

' 其中T 为系统的动能，包括小球的转动的动能，导轨转动的动能等，V 为系统的势能，包括重力势能弹性势能等等，能量耗散函数为R ，()

T

k q q q q (21)

为广义坐标向量，其中k 代表系统的自由度，即完全描述系统运动特性需要的坐标数目，关于自由度在下文会具体分析，u 为作用于系统的外力。

以下为各个变量所表示的物理意义，M ：导轨的质量，g ：重力加速度 r ：小球的半径 b I ：球的惯性力矩，w I ：杆的惯性力矩，x ：球的相对横坐标，y ：球的相对纵坐标，ψ：小球相对于导轨的转角，a ：导轨与水平线的夹角，球杆系统受力分析如下：

图1. 球杆受力分析

由图可知，小球的广义坐标向量有3个，即是小球坐标x ，小球转动角ψ，杆的转角α，而小球的速度和导轨转动的角速度可以用?

x 和?

a 来表达，因此，通过这三个变量，就能完全描述系统的各个状态，但是在实际中，由于小球是沿着导轨只做滚动不做滑动的，因此就有ψR x =，也就是说小球的转动角ψ和小球位移是一一对应的，那么转动角ψ就可以和x 合并到一起分析，这样球杆系统有两个广义坐标x q =1，a q =2。

要想计算系统的动能，就必须知道小球和导轨的转动惯量，根据资料得小球和导轨的转动惯量如下：

2521032.45

2

kgm mR I b -?==

21402.0kgm I w =

导轨是一个重心在重心的均匀旋转体，其角速度为?

a ，因此，导轨的动能为:

2

2

1?=a I T w w

而小球的动能包括两部分，第一部分是小球沿着导轨运动的动能，用2

2

1

mv

来计算，第二部分是小球的旋转动能，用小球的转动惯量和角速度来计算，公式如下：

22

2

121mv I T b b b +=ω

其中，小球的转动角速度b ω要考虑到导轨的运动，因此要加上导轨的角速度?

a

??

+=a R

x b ω

222

2

2??

??++=a R

a x R x b

ω 因为导轨转动角速度?

a 速度很小，因此忽略后两项得：

2

2

2

R

x b

?=ω

而式中v 是小球的线速度，根据资料可以用向量积来表示：

r w v v ?+='

其中'v 是小球相对于导轨的线速度，其数值等于?

-x ，负号是指方向与规定的正方向相反，ω指的是导轨的角速度，即?

a ，r 是小球的质心在坐标系中的位置向量，计算式如下：

?????

???

????---=??????????-???????????+????????????-=?????000

000a x R a x R x a x v 由于小球的半径R 很小，同时导轨转动的角速度?

a 也很小，因此上式可以近似等于：

?????

?

??????--≈??0a x x v

2

222??+=a

x x v

将以上结果带入小球动能计算公式中，得小球的动能为：

????

???????????

??+???? ??+=???222221R x I a x x m T b b 而前文已经得出导轨的动能w T ，因此两式相加得系统的动能T 为：

????

??????+????? ??+???? ??+=+=????2222

221a I R x I a x x m T T T w b w b 而系统的势能为小球的重力势能，以导轨水平位置为零点，得势能V 为：

a mgx V sin -=

以上，系统的动能和势能均已求出，下面带入拉格朗日方程进行求解，由于篇幅原因，推导过程在本文中不再赘述，带入拉格朗日方程后得到关于导轨倾角

α的拉格朗日方程为：

()τ=-++?

???a mgx a x mx a I mx w

cos 22

小球坐标x 的运动方程为：

a mg a mx x R I m

b sin 22=-???

?

?+???

球杆系统的运动方程为：

a ma mxa x R I m

b sin '2

2=-??? ?

?+??

()τ=-++?

???a mgx a x mx a I mx w

cos 22

2.3拉格朗日模型线性化及状态空间表达式求取

以上为通过拉格朗日方法建立起来的系统运动方程，但是由于运动方程存在非线性环节，无法直接转化为状态空间表达式，需要通过在稳定工作点附近做近似线性化处理来消除，根据实际情况来说，我们一般考虑的是小球在一个确定点附近的稳定性，假设这个点为0x ，那么这个点就是系统的静态工作点，动能计算

式中2m x 项可以近似认为是2

0mx ，因此2m x 就可以近似认为是常数，运动方程中的相关倒数项就消失了，因此运动方程近似为：

a mg x R I m

b sin 2=??? ?

?

+??

()τ=-+?

?a mgx a I mx w

cos 2

进一步进行简化，由于导轨倾角a 一般较小，因此近似取a a ≈sin ，1cos ≈a ，

得到最终的运动方程为：

mga x R I m b =??? ?

?

+??2

()

τ

=-+?

?mgx a I mx

w 2

根据资料，小球半径m R 01899.0=、小球质量kg m 03.0=

选取状态变量x x =1，?

=x x 2，a x =3，?

=a x 4，带入参数即得到以下状态空

间表达式：

F

X X X X X X X X ?

????

?

??????+????????????????????????=???????

?

????????495.3000000873

.181000070000104321.4.3.

2.1

????

?

????????????

?=??????43212101000001X X X X Y Y

3.系统稳定性分析

系统的稳定性是系统的根本特性，一个不稳定的系统几乎是无法使用的，为了分析一个系统，首先应该分析其稳定性，根据常识，球杆系统在不给予控制时是不稳定的，因为当导轨的倾角是固定时，小球将会在重力的作用下一直运动到

导轨底部，即使导轨停留在水平位置时，小球可能处于静止状态，但是当导轨受到一个微小的扰动后，小球就无法回到原来的位置了。如果导轨既不是固定倾角，又不是水平状态，而是无规律摆动的话，那么小球更不可能稳定在一个位置，要使小球能够停在任意位置，必须对导轨的倾角施加控制，下面通过Matlab来分析系统在有初始条件下的稳定性。

3.1有初始状态下求取系统响应曲线

设置小球初始位移为0.5，导轨初始倾角为5°，进行仿真，程序如下：

A=[0 1 0 0;0 0 7 0;0 0 0 1;18.873 0 0 0];

B=[0 0 0 3.495]';

C=[1 0 0 0;0 0 1 0];

D=0;

T=0:0.005:10;

U=zeros(size(T));

X0=[0.5 0 5 0];

lsim(A,B,C,D,U,T,X0);

响应曲线如下：

图3.有初始状态响应曲线

由上图可以看出，小球的位移和导轨的倾角在9秒钟以后都开始发散，从实际情况来说，因为小球有初始位置，导轨也有初始倾角，那么小球必然在重力的作用下加速下滚，而导轨的倾角也将越来越大，在一段时间后，小球的位置和导轨的倾角都将发散，这也从侧面说明所建模型是合理的。

3.3稳定性判断并求取零极点分布图

以上的分析是根据常识分析的，下面，通过严格的判断方法来分析系统的稳定性，系统稳定性分为外部稳定和内部稳定，外部稳定是指输入和输出之间的关系，由系统传递函数分母的特征根决定，但是当系统不稳定的极点对消时，有可能出现假稳定的情况，而内部稳定是指系统内部变量能够趋于稳定，而由内部稳定可以推出外部稳定，反之则不行，因此现在只判断内部稳定，内部稳定的充要条件是特征矩阵的所有特征根具有负的实部，这样系统的能量就逐渐减小，最终达到稳定条件，如果特征值有正的实部，那么能量就不能衰减，导致系统不稳定。

求系统特征根分布的程序如下：

A=[0 1 0 0;0 0 7 0;0 0 0 1;18.873 0 0 0];

B=[0 0 0 3.495]';

C=[1 0 0 0;0 0 1 0];

D=0;

eig(A)

pzmap(A,B,C,D);

求得的特征根为

ans =

-3.3903 + 0.0000i

0.0000 + 3.3903i

0.0000 - 3.3903i

3.3903 + 0.0000i

零极点坐标图如下：

图4.系统零极点图

由系统零极点坐标图可知，系统有一个极点位于系统右半平面，并且有一对虚轴上的共轭虚根，因此，系统是不稳定的。

4.系统能控性判别

由上文可知，系统是不稳定的，我们知道，稳定是一个系统能够使用的前提，一个不稳定的系统是没有意义的，但是在实际使用过程中，我们可以通过状态反馈或者输出反馈等方法来重新配置系统的极点，使系统达到稳定，而重新配置极点的充要条件是：一、系统完全能控，这样所有的特征根都可以重新配置。二、系统部分可控，但是不可控极点拥有负的实部，而拥有正实部的极点是可控的，这种情况下我们依然可以通过配置可控极点使系统稳定。

一般来说，判断能控性的方法主要有两种，一种是代数判据，一种是模态判据，下面分别通过两种方法判断系统能空性。

4.1代数判据

对于一个线性定常系统，系统完全能控的充分必要条件是能控判别矩阵

n n n c B A AB A U ?-=]....[1满秩，而Matlab 中可以通过ctrb 命令轻松求得能控

矩阵，配合求矩阵秩的命令rank 可以判断矩阵是否能控。

能控性判别程序如下：

A=[0 1 0 0;0 0 7 0;0 0 0 1;18.873 0 0 0]; B=[0 0 0 3.495]'; C=[1 0 0 0;0 0 1 0]; D=0;

Uc=ctrb(A,B); N=rank(Uc) if N==4

disp('系统完全能控'); elseif N<4

disp('系统不完全能控') end

结果为

N = 4 系统完全能控

由此可知，系统是完全能控的，也就是说我们可以通过状态反馈随意配置系统的所有极点。

4.2模态判据

模态判据分为特征根两两互异和特征根有重根两种情况，当特征根两两互异

时，模态判据表示为，把系统化为对角标准型：

u b b b b x x p

n np n p n n n ??????????

???????+?

?????

???

???=...........................

111121λλλ

以后，B 矩阵，没有全为零的行即表示系统完全可控。

当特征根有重根时，模态矩阵表示为，把系统化为约当标准型，这里假设系统有p 重根λ和p n -个互不相等的特征根，约当标准型如下：

Bu x B A x +???

??

?=?

其中A 矩阵为p p ??

???????????λλλ1......1 B 矩阵为)

()(21...

p n p n p n -?--????

???

???????λλλ 化为约当标准型以后，只要最后一行不全为0，那么矩阵完全可控。

根据前文可知，本系统的四个特征根两两互异，因此，只需化成对角标准型即可，Matlab 程序如下：

A=[0 1 0 0;0 0 7 0;0 0 0 1;18.873 0 0 0]; B=[0 0 0 3.495]'; C=[1 0 0 0;0 0 1 0]; D=0;

sys=ss(A,B,C,D);

sysc=canon(sys,'modal') 运行后得结果：sysc =

a =

x1 x2 x3 x4 x1 -3.39 0 0 0 x2 0 4.996e-16 3.39 0 x3 0 -3.39 4.996e-16 0 x4 0 0 0 3.39

b =

u1 x1 0.4578

x2 0.6808 x3 -0.0925 x4 -0.4578

c =

x1 x2 x3 x4

y1 -0.3429 0.06151 0.4527 -0.3429 y2 -0.563 -0.101 -0.7434 -0.563

d = u1 y1 0 y2 0

即对角标准型为：

u e e x ????????????--+????????????----=?4578.00925.06808.04578.039.3000016996.439.30039.316996.4000039

.3 x y ??

????------=563.07434.0101.0563.03429.04527.006151.03429.0 可以看出，B 矩阵没有任意一行完全为零，因此系统能控。

4.3可控性与可稳定性

前文曾经提过，如果系统所有的极点都是可控的话，那么系统可以通过状态观测器任意配置极点，从本文来看，即可以达到使系统由不稳定配置到稳定，那么是不是不完全能控的系统就不能配置到稳定呢？答案是否定的，事实上，在部分可控的系统中，如果所有不可控的极点，只要其满足()0Re

5.系统极点配置

由上文可知，系统是完全可控的，即所有极点都可以通过状态反馈配置到我们想要的位置，以获得期望的动态性能，下面从极点配置方法入手，对极点配置部分进行研究。

5.1极点配置方法

一般来说，我们对于系统的极点配置方法主要有三种，分别是状态反馈，输出反馈和用PID 配置极点，这三种方法从理论上来说都可以对极点进行配置，下边对三种方法进行分析。 5.1.1状态反馈原理

图5.状态反馈原理图

根据原理图，状态反馈是把系统的状态通过一个转换矩阵K 引到输入端，这样输入就变成了kx r u -=，带入后即有()Br x BK A x +-=?

这样，系统的特征矩阵即变成了()Bk A -，而系统的特征矩阵决定了系统的特性，这样我们就通过配置反馈矩阵k 来将()Bk A -的闭环极点设置到我们想要的位置。

5.1.2输出反馈原理

图6.输出反馈原理图

跟状态反馈类似，输出反馈是将输出通过一个转换矩阵F 引到输入端，这样输入就变成了，Fy r u -=，带入后有()Br x BFC A x +-=?

，系统特征矩阵变成了

()BFC A -，同样我们可以通过配置矩阵F 来配置系统的极点。

5.1.3PID 配置极点原理

图7.PID 配置极点

PID 配置极点的方法是将PID 传递函数写成

S

S K K S K d i p 2

++的形式，带入闭

环中得到闭环传递函数，这时把闭环传递函数的分母等于()()()321λλλ---S S S ，

321,,λλλ则是期望的极点，用待定系数法即可求得PID 参数。

5.1.4三种反馈对比

这三种方法都可以配置极点，那么如何选择呢，事实上，输出反馈不如状态反馈适应能力强，因为一般系统，输出的数量要少于状态的数量，比如本文中的系统，就是4个状态，2个输出，这就代表，用状态反馈有4个自由度，而用输出反馈仅有两个自由度，由于输出反馈自由度要小，能用状态反馈配置的时候，输出反馈并不一定能够配置成功，而PID 配置极点的方法由于结构限制，对于二阶一下系统较为实用，而高阶系统则不好用，因此，本文选用状态反馈进行极点配置。

5.2.用状态反馈进行极点配置

闭环极点的选择首先要保证系统的稳定性，即期望的特征根应有负的实部，其次考虑到系统的动态特性，一般选择的极点是原极点的两到三倍，因此我选择极点如下：

61-=s 72-=s J s 664,3±-=

Matlab 程序如下：

A=[0 1 0 0;0 0 7 0;0 0 0 1;18.873 0 0 0]; B=[0 0 0 3.495]'; C=[1 0 0 0;0 0 1 0]; D=0;

G=ss(A,B,C,D); figure(1) step(G)

M=[-6 -7 -6+6j -6-6j];

K=acker(A,B,M)

G1=ss(A-B*K,B,C,D);

G2=tf(G1);

H=eig(A-B*K)

figure(2)

step(G2)

得到控制器参数为K = 129.0052 58.8596 77.2532 7.1531 下面，对比配置前后的响应曲线，配置以前的响应曲线如下：

图8.极点配置前阶跃响应图

极点配置以后的响应曲线如下：

图9.极点配置后阶跃响应图

由极点配置以后的响应曲线看出，小球位移和导轨的倾角在1.2S以后均趋于稳定，这表示小球可以控制到导轨上的任意位置，极点配置是成功的。

为了进一步验证极点配置的有效性，我对极点配置以后的系统在有位移为0.5，倾角为5°的情况下进行了仿真，仿真程序如下：

A=[0 1 0 0;0 0 7 0;0 0 0 1;18.873 0 0 0];

B=[0 0 0 3.495]';

C=[1 0 0 0;0 0 1 0];

D=0;

K =[129.0052 58.8596 77.2532 7.1531];

T=0:0.005:10;

U=zeros(size(T));

X0=[0.5 0 5 0];

lsim(A-B*K,B,C,D,U,T,X0);

得到的曲线如下：

图10，极点配置后有初始状态的响应

在配置极点以前，同样在0.5米5°的初始条件下，仿真曲线如下：

图11.极点配置前有初始状态的响应

根据这两个曲线的对比来看，极点配置很好的解决了在有初始条件下系统的稳定性。

6.可观性分析及带状态反馈的状态观测器的设计

6.1能观性分析

本文系统有两个输出，分别是位移X 和倾角α，这两个状态可以直接从输出中获得，而系统的状态变量有4个，分别是X ，?

X ，α和?α而?X 和?

α是不能从输出中直接获得的，而我们在实际使用中往往是需要这些状态的，这就需要设置一个全维状态观测器来观测无法直接获得的变量，而全维状态观测器能够设置的充分条件是系统完全能观，下面我讨论系统的能观性，类似于能控性判别，我们仍然使用代数判据和模态判据两种方法来判断能观性。 6.1.1代数判据 Matlab 程序如下：

A=[0 1 0 0;0 0 7 0;0 0 0 1;18.873 0 0 0]; B=[0 0 0 3.495]'; C=[1 0 0 0;0 0 1 0]; D=0;

Uo=obsv(A,C); N=rank(Uo)

if N==4

disp('系统完全能观') elseif N<4

disp('系统不完全能观')

end

运行程序得结果为N=4 系统完全能观

这个结果说明系统是完全能观的，那么我们就能通过状态观测器来得到无法从输出得到的那两个变量。 6.1.2模态判据

能观性模态判据与能控性相似，只不过判断标准由B 矩阵没有任何一行全为零变为了C 矩阵没有任何一列全为0，有前文可轻松判断出模态判据下，系统也是完全能观的，在此不再赘述。

6.2能观性与可检测性

类似于能控性与能稳定性的条件，一个不完全能观的系统，只要其不能观极点满足()0Re

那必然是能检测的，因此不对能检测性做判断。

6.3全维观测器原理

为了能够通过状态反馈配置系统极点，就必须知道系统中各个状态变量的变化曲线，而实际系统中往往不是每个状态都能直接测得的，以本系统为例，X1和X3分别是小球的位移和导轨的转角，这两个量都是可以作为输出直接测量的，而X2和X4分别是小球的速度和导轨的角速度，这两个量是无法直接测量的，因此要用状态变量进行状态反馈，必须构建全维状态观测器，并从观测器中提取各个状态变量，因此，状态反馈和状态观测器是相互联系的。

状态观测器的原理即状态重构，即在实际系统结构已知的条件下，通过建模手段在计算机中构建一个机理模型，并且给这个模型输入与实际系统相同的输入信号，这样，就相当于在计算机中模拟了整个实际系统，从模拟的过程中我们就可以得到系统各个状态的观测值，并利用这些状态通过状态反馈矩阵K来重新配置极点。

6.4全维状态观测器结构

带状态反馈的全维状态观测器的结构如下：

图12.带状态反馈的全维状态观测器原理图

状态观测器上半部分是实际的系统，下半部分是在计算机中进行仿真的观测器，因此用虚线方框框住，代表这部分是虚拟的，因为仿真总会有误差，而且实

际系统容易受到干扰，为了避免误差逐步累加，因此用实际输出和观测器输出的差值来校正观测器。

根据原理图可得，虚线方框中状态观测器的状态空间表达式为：

??

?

??-++=?

^

^^

X c y L BU X A X 即()Ly Bu X LC A X ++-=∧

?∧

6.5全维状态观测器设计

根据上文，我们只要将（A-LC ）的极点配置到想要的位置即可求得状态观测器的L 矩阵，而观测器的基本要求是反应速度快，这样才能更快的跟踪原系统，因此一半选择观测器极点为原极点的4-6倍，现设置极点为：

121-=p ，112-=p ，j p 6114,3±-=

状态观测器求取的Matlab 程序如下：

A=[0 1 0 0;0 0 7 0;0 0 0 1;18.873 0 0 0]; B=[0 0 0 3.495]'; C=[1 0 0 0;0 0 1 0]; D=0;

p=[-12 -11 -11+6j -11-6j]; L=place(A',C',p)

求得L =

23.5469 136.8709 4.3796 70.3498 -8.0089 -82.1087 21.4531 117.8992 而状态反馈矩阵在前文已经求得，即 K =

129.0052 58.8596 77.2532 7.1531

6.6全维状态观测器Simulink 仿真

前文已经求得带状态反馈的全维状态观测器的全部参数，在Simulink 中搭建出模型并进行仿真，由于系统为4阶系统，并且有两个输出，结构较为复杂，因此K 矩阵和L 矩阵均使用了Subsystem 来简化，模型截图如下：