第七章matlab.ppt

1、第7章 控制系统的状态空间设计法,状态空间设计法是一类以系统状态空间模型（即状态空间表示式）为基础的设计方法；其任务在于设计控制器，寻求改善系统性能的各种控制规律，以保证系统的各项性能指标要求都得到满足。MATLAB环境最适于矩阵运算，这使得进行状态空间设计非常方便，本节将详细讨论控制系统的状态空间设计问题并给出MATLAB实现。,7.1 极点配置,线性系统的动态性能，如系统稳定性，时域分析中的超调量、过渡时间等指标，主要取决于系统的极点位置。因此，作为系统性能指标的一种形式，往往是给定一组期望极点，或者根据时域指标转换成一组等价的期望极点。极点配置问题，就是通过选择反馈增益矩阵，将闭环系统的

2、极点恰好配置在根平面上所期望的位置，以获得所希望的动态性能。,7.1.1 单输入系统的极点配置,对于SISO系统来说，极点配置问题描述如下： 设控制系统的状态方程模型为 (7-1) 当引入状态反馈后，系统的控制信号为，这里是系统外部的输入；是行向量，称为反馈矩阵；此时闭环系统状态方程模型为： (7-2),可以证明，当系统状态完全可控时，则可以通过状态反馈将系统的极点配置到复平面的任何位置。 极点配置有两种方法：第一种方法是采用变换矩阵，使系统具有期望的极点，从而求出矩阵；第二种方法基于Caylay-Hamilton理论，通过矩阵特征多项式，可求出（这种方法称为Ackermann公式）。,1、矩

3、阵变换法 第1步：系统通过状态反馈成立 (7-3) (7-4) 式中，为期望的闭环极点（实数极点或共轭复数极点）。,若完全能控，必存在非奇异变换 (7-5) 式中 (7-6) 其中， 为可控性矩阵,阵定义为,能将化成能控标准I型 (7-7) 式中,加入状态反馈增益阵 (7-8) 可求得对的闭环状态方程表达式 (7-9) 式中,闭环特征多项式为 （7-10） 使闭环极点与给定的期望极点相符，必须满足 由等式两边同次幂系数对应相等，可解出反馈阵各系数 (7-11) 于是得,最后，把对应于 的 ，通过如下变换，得到对应于状态 的K。 (7-12) 这是由于 的缘故。,以上为极点配置的理论基础，因此，

4、用MATLAB进行极点配置的步骤为: 第1步 检查系统的可控性，当 时系统可控；在MATLAB中可用rank(ctrb(A,b)实现； 第2步 确定系统矩阵A的特征多项式系数 在MATLAB中，可用poly函数实现； 第3步 确定变换矩阵T cam为能控性判别矩阵，在MATLAB中可用ctrb(a,b)实 现，w在MATLAB中可用hankel()函数来实现。,第4步 确定期望特征多项式系数 第5步 求增益矩阵,2、Ackermann公式 由上面讨论知，闭环系统 特征方程为 由此定义 则对系统矩阵 有 根据Cayley-Hamilton理论， 成立，即,将 代入，可得 在MATLAB中，可很方

5、便地利用上述两种方法进行极点配置设计，但利用MATLAB控制系统工具箱提供的place和acker函数进行极点配置设计，可免去繁琐的过程，因此一般可直接采用这两个函数进行设计。,功能： SISO极点配置 格式： K = acker(a,b,p) 说明： acker函数利用Ackermann公式计算反馈增益矩 阵K，使采用全反馈的单输入系统具有指定的闭环极 点，即,3、acker函数,4、 place函数,4、 place函数 功能：MIMO极点配置 格式：K = place(A,B,p) K,prec,message = place(A,B,p) 说明：place函数利用Ackermann公式

6、计算反馈 增益矩阵K，使采用全反馈的多输入多输出系统具有指定的闭环极点 ，即 K,prec,message = place(A,B,p)同时返回系统闭环极点与希望极点p的接近程度prec。prec中的每个量的值为匹配的位数。如果系统闭环极点的实际极点偏离希望极点10%以上，则message将给出警告信息。,系统可任意极点配置的前提条件是系统必须可控，在进行极点配置前必须判断系统的可控性，判别系统可控性程序如下 function m = controllble(a,b) nc = rank(ctrb(a,b); %求可控矩阵的秩，ctrb(a,b)为系统的可控矩阵 n = length(a);

7、if n =nc %判断系统可控矩阵是否满秩 flagc = system is completely state controllable m = 1; else flagc = system is not completely state controllable m = 0; end,例7-1 已知系统的传递函数为 希望极点为 ，试设计状态反馈矩阵 。 下面由命令行分别作两种方法的极点配置计算，即矩阵变换法与Ackermann极点配置算法。,7.1.2 极点配置实例分析,1、矩阵变换法 程序 num = 1; den = 1 18 72 0; a,b,c,d = tf2ss(num,de

8、n); m = controllble(a,b); %判别系统的能控性。 if m = 1 %系统满秩，系统状态是完全能控的，可以通过状态反馈极点配置。 fy = poly(a) fq = conv(1 10,conv(1 2+2j,1 2-2j) w = hankel(fy(length(fy)-1:-1:2);1) cam = ctrb(a,b) T = cam*w i = length(fy):-1:2 diffa = -(fy(i)-fq(i) K = diffa*inv(T) else message(This system not controllable ,can not pol

9、e allocation) end,2、Ackermann法 程序 num = 1; den = 1 18 72 0; a,b,c,d = tf2ss(num,den) m = controllble(a,b) if m = 1 %系统满秩，系统状态是完全能控的，可以通过 状反馈极点配置。 fy = poly(a) %计算原系统特征多项式 fq = poly(-10,-2+2i,-2-2i) %计算配置极点的特征多项式 thta = polyvalm(fq,a) k = 0 0 1*inv(ctrb(a,b)*thta else message(This system not controll

10、able ,can not pole allocation) end,MATLAB控制系统工具箱提供的place和acker函数，可直接采用这两个函数进行极点配置设计，免去了繁琐的编程计算过程。MATLAB程序如下： num = 1; den = 1 18 72 0; a,b,c,d = tf2ss(num,den) %判别系统的能控性。 m = controllble(a,b) if m = 1 %系统满秩，系统状态是完全能控的，可以通过状态反馈极点配置。 p = -10 -2+2j -2-2j; K = acker(a,b,p) %或 K = place(a,b,p) else messa

11、ge(This system not controllable ,can not pole allocation) end,系统运行结果为： k = -4.0000 -24.0000 80.0000 系统的特征值为 eig(a-b*k) ans = -10.0000 -2.0000 + 2.0000i -2.0000 - 2.0000i 例7-2 设系统的传递函数为 试设计状态反馈控制器，使闭环系统的极点 为-2， 10，-10。,解： 下面由命令行分别作两种方法的极点配置计算，即非奇异变换法与Ackermann极点配置算法。 程序 num = 10; den = conv(1 0,conv(

12、1 1,1 2); a,b,c,d = tf2ss(num,den) %判别系统的能控性。 m = controllble(a,b) if m = 1 %系统满秩，系统状态是完全能控的，可以通过状态反馈极点配置。 p = -2 -10 -10; k = place(a,b,p) else message(This system not controllable ,can not pole allocation) end,系统运行结果为： ? Error using = place Cant place poles with multiplicity greater than rank(B).

13、观测器极点为重根时，提示出错信息。由此可见place()命令不适于进行极点为重根的极点配置，可另取近似重根计算，这样就可避免上述错误出现。对于有重根的情况可直接用acker()命令计算 K = acker(a, b, p) k = 1.0000 4.0000 4.0000 系统的特征值为 eig(a-b*k) ans = -2.0000 -1.0000 + 1.0000i -1.0000 - 1.0000i 无上述重根问题出现，但是acker()命令只适用于SISO系统，应用时要注意。,Ackermann法: num = 10; den = conv(1 0,conv(1 1,1 2); a,

14、b,c,d = tf2ss(num,den) m = controllble(a,b) if m = 1 %系统满秩，系统状态是完全能控的，可以 通过状态反馈极点配置。 fy = poly(a) %计算原系统特征多项式 fq = poly(-10,-2+2i,-2-2i) 计算配置极点的特征多项式 hemil = polyvalm(fq,a) k = 0 0 1*inv(ctrb(a,b)*hemil else message(This system not controllable ,can not pole allocation) end,状态观测器，其原理结构如图7-1所示。 当观测器的

15、状态 与系统实际状态 不相等时，反映到它们的输出 与 也不相等，于是产生一误差信号 ，经反馈矩阵馈 送到观测器中每个积分器的输入端，参与调整观测器状态 ，使其以一定的精度和速度趋近于系统的真实状态 。,7.2 观测器设计及实例分析,根据图7-1a可得状态观测器方程 即 (7-13) 式中，为状态观测器的状态矢量，是状态 的估值； 为状态观测器的输出矢量； 为状态观测器的输出误差反馈矩阵。 根据式(7-13)，可将状态观测器表示成图7-1b。从图中看出，它有两个输入，一个是待观测系统的控制作用 ，一个是待观测系统的输出 。它的一个输出就是状态估值 。,系统可设计状态观测器的前提条件是系统必须可观

16、，在设计状态观测器前必须判断系统的可观性，判断系统可观性程序如下: function m = observable(a,c) %判断系统的可观性 no = rank(obsv(a,c); %判断系统可观矩阵的秩 n = size(a,1); if n= no %判断系统可观矩阵是否满秩 flago = system is complete observable m = 1; else flago = system is not complete observable m = 0 end,例7-3 线性系统 给定观测器极点为 ，计算观测器反馈矩阵,程序 a = 0 1;-2 -3; b = 0;

17、1; c = 2 0; d = 0 m = observable (a,c) if m = 1 %系统满秩，系统状态是完全能观的，以设计状态观测器。 p = -10 -10; k = acker(a,c,p) else message(This system not observable ,can not design observer) endG=acker(a,c,-10 -10) MATLAB程序运行结果为 G = 8.5000 23.5000,Ackermann法 a = 0 1;-2 -3; b = 0;1; c = 2 0; d = 0 m = observable (a,c) i

18、f m = 1 %系统满秩，系统状态是完全能观的，可以设计状态观测器。 fy = poly(a) %计算原系统特征多项式 fq = poly(-10 -10) %计算配置极点的特征多项式 pthi = polyvalm(fq,a) k = 0 0 1*inv(ctrb(a,b)*pthi else message(This system not controllable ,can not pole allocation) end,由观测器的状态方程可以看到，观测器的系数矩阵的 计算是状态反馈矩阵 的转置计算，即观测器方程的对偶系统的状态反馈，因此MATLAB函数place()和acker()同

19、样可以用于状态观测器反馈矩阵 的计算。,例7-4 已知系统 设计状态观测器使其极点为 。,程序 a = 1 0;0 0; b = 1;1; c = 2 -1; d = 0 m = observable(a,c) if m = 1 %系统满秩，系统状态是完全能观的，可以设计状态观测器。 p = -10 -10; G = acker(a,c,p) else message(This system not obvserble ,can not design observer) end,程序运行结果为 G = 60.5000 100.0000 观测器方程为 或者,模拟结构图如图7-2所示。,状态观测器

20、解决了受控系统的状态重构问题，可使状态反馈系统得以实现。 图7-3是一个带有全维状态观测器的状态反馈系统 图7-3 带状态观测器的状态反馈系统。,7.3 基于观测器的状态反馈控制器设计及实例分析,设能控能观的受控系统 为 (7-14) 状态观测器 为 (7-15) 反馈控制律为 (7-16),将式(7-16)代入式(7-14)和式(7-15)整理或直接由结构图得整个闭环系统的状态空间表达式为 (7-17) 写成矩阵形式为 ， 即 (7-18),闭环系统的极点包括 直接状态反馈系统 的极点和观测器 的极点两部分。但二者独立，相互分离。 设状态估计误差为 ，引入等效变换： (7-19) 令变换阵为

21、 （7-20）,经线性变换后的系统 为,或者展开成 其等效结构如图7-4所示。,由于线性变换不改变系统的极点，因此， 有 (7-23),例7-5 设受控系统的传递函数为 用状态反馈将闭环系统极点配置为 设计实现上述反馈的全维观测器，设其极点为(-10，-10)。,程序 num = 1; den = conv(1 0,1 6); a,b,c,d = tf2ss(num,den) %判别系统的能控性。 m = controllble(a,b); %判别系统的能控性。 n = observable(a,c) if m =1 K = acker(a,b,p1) G = acker(a,c,p2) el

22、se message(This system not controllable ,can not pole allocation) end,程序运行结果为： 全维观测器方程为 闭环系统结构如图7-5所示。,例7-6 单倒立摆系统原理图如图3-27所示。长度为 、质量为 的单倒立摆，用铰链安装在质量为 的小车上，小车受执行电机操纵，在水平方向施加控制力 ，相对参考系产生位移 。为简化问题并保留实质，忽略执行电机惯性，摆轴、轮轴、轮与接触面之间的摩擦及风力。若不给小车施加控制力，倒立摆会向左或向右倾倒，故是一个不稳定系统。目的是当倒立摆出现偏角 以后还能保持在垂直位置。要求建立该系统数学模型，用状

23、态反馈配置系统极点，运用全维状态观测器实现状态反馈。,解 为了分析与设计的需要，首先建立倒立摆的运动方程并将其线性化，注意在此处设摆杆长度为 ，而图3-27所示中的摆杆长度为 。设小车瞬时位置为 ，摆心瞬时位置为（ ）。在 作用下，小车及摆均产生加速运动，根据牛顿第二定律，在水平直线运动方向的惯性力应与 平衡。于是有： 即 绕摆轴旋转运动的惯性力矩应与重力矩平衡，于是有 即,由于控制目的是保持倒立摆直立，在施加适宜的条件下，假定 、 均接近于零是合理的， 此时 ， ，且可忽略项 ，于是有 联立求解可得 消元可得四阶系统微分方程,选取小车的位移 及其速度 、摆的角位移 及其角速度 作为状态变量，

24、 作为输出变量，考虑恒等式 ， ，可列出系统的状态空间表达式为 式中,假定系统参数为： ， ， ， ，则MATLAB程序如下： M = 1; m = 0.1; l = 1; g = 9.8 A = 0 1 0 0;0 0 -(m*g/M) 0;0 0 0 1;0 0 (M+m)*g/(M*l) 0; b = 0;1/M;0;-(1/M*l); c = 1 0 0 0; m = controllble(A,b) n = observable(A,c) if m =1 0 -2; B = 0;1; C = 1 0; D = 0; Q = 1 0;0 2; R = eye(1); p = care(

25、A,B,Q,R) 求得的本系统的代数Riccati方程的解为 p = 2.8284 1.0000 1.0000 0.8284,五、降维观测器,（2）、求解系统的最优控制信号 根据最优控制信号的表达式： 得到系统的代数Riccati方程的解之后可以求出最优控制信号的变化曲线，其中，状态反馈矩阵K等于 ，紧接上一步，输入以下代码：,K = -inv(R)*B*p 求得K为 K= -1.0000 -0.82842 T = 0:0.005:15 %仿真时间 sys = ss(A+B*K,B,C,D); %系统状态方程描述 y,T,xt = step(sys); n = length(T); %求解最优

26、控制信号u(t) m = length(R); for i = 1:m for j = 1:n u(j,:,i) = K(i,:)*(xt(j,:,i); end end plot(T,u(:,:,1) %图形绘制 title(Linear Quadratic Controller Output u(t) xlabel(Time(sec) ylabel(Value); grid;,得到系统的最优控制信号如图7-7所示。,（3）、系统对输入的响应曲线 输入代码如下： step(A+B*K,B,C,D,1,T) 这里求解系统的闭环参数矩阵时用的是A+B*K而不是A-B*K，这是因为在计算状态反馈矩

27、阵K时已经包含负号。响应结果如图7-8所示 本例也可直接利用k,p,e = lqr(A,B,Q,R)直接获得最优控制状态反馈矩阵K。,例7-8 某倒立单摆系统如图7-9所示，其状态方程为：,其中小车的质量为 ，倒立单摆的质量 为，小车的摩擦系数 为，端点与倒立单摆质心的距离为 ，倒立单摆的惯量为 ，输入量 是施加在小车上的外力，四个状态变量分别是小车的坐标 、 ，倒立单摆的垂直角度 、 。输出的被控量分别是小车的坐标 和倒立单摆的垂直角度 。试根据误差指标J最优意义下最优的规则设计线性二次型最优控制器和相关的参考输入以及观测器，满足以下指标 （1）、输出量 和 的过渡过程时间小于2s。 （2）

28、、输出量的上升时间小于0.5s （3）、输出量 的超调显小于 ( ）,求解过程： 本例的解题步骤分为以下几步： （1）、分析原系统的开环阶跃响应 首先要求得其开环系统的特征值，判断其稳定性。然后根据其阶跃响应曲线分析当前的运动情况与期望性能指标之间的差距，确定校正的手段。在MATLAB中输入以下代码： %二次型最优控制例子 %清除所有内存变量 clear all %系统参数初始化 M = 0.5; m = 0.2; b = 0.1; i = 0.006; g = 9.8; l = 0.3;,%系统状态方程描述 p=i*(M+m)+M*m*l2; A = 0 1 0 0;0 -(i+m*l2)*

29、b/p (m2*g*l2)/p 0;0 0 0 1; 0 -(m*l*b)/p m*g*l*(M+m)/p 0; B = 0;(i+m*l2)/p;0;m*l/p; C = 1 0 0 0; 0 0 1 0; D = 0;0; %求解系统的特征值 p = eig(A) t = 0:0.01:1; step(A,B,C,D,1,t) 求得系统的特征值为 p = 0 5.5651 -0.1428 -5.6041,系统有一个右半平面的极点，因此不稳定，必须加入校正装置。此时系统的阶跃响应曲线如图7-10所示。,（2）、线性二次型最优控制器设计 设计线性二次型最优控制器关键是选择加权矩阵Q。一般说来，

30、Q选择的越大，系统达到稳态所需的时间越短，当然，还要实际的系统允许。 首先选择 ， ， 然后根据实际情况进行调节。输入以下代码： x = 1; y = 1; Q = x 0 0 0; 0 0 0 0;0 0 y 0; 0 0 0 0 ; R = 1; %求解线性二次型最忧状态反 K = lqr(A,B,Q,R) %求解系统闭环状态力程 Ac = (A-B*K); Bc = B; Cc = C; Dc = D; %辅出系统阶跃仿真,T = 0:0.02:5; U = ones(size(T); Y,X = lsim(Ac,Bc,Cc,Dc,U,T); plot(T,Y(:,1),.-,T,Y(:,2); title(Inverted Pendulum LQ Step Response); xlabel(Time-sec); ylabel(Respo