9 最优控制
用经典控制理论设计系统时，是根据给定的频域和时域指标，通过选择适当结构和参数的校正装置进行调节，本质上是一个试凑的方法，设计质量的好坏很大程度上依赖于设计人员的实践经验。而最优控制的设计目标是要选择适当的控制的规律，使控制系统在严格的数学基础上实现系统的性能品质在某种意义下是最优的。最优控制是现代控制理论的一个重要组成部分。
本章主要讨论最优控制问题的基本概念。二次型最优控制问题以及最优观测器设计。
9．1最优控制问题实例
9．1．1最优控制问题实例
设有一盛放液体的连续搅拌槽，如图9.3所示。槽内装有搅拌器使槽内液体经常处于完全混合状态。若槽内液体初始温度为0℃，现需将其温度经1小时后升高到40℃。为此在
入口处送进一定量的温度为u(t)的加热液体。 （为保持槽内液体恒定，出口处的流出量等于 入口处的流入量），试寻找加热液体温度u(t)的变化规律，使槽中的液体温度经1小时后上升到40℃，且散失热量最小。
因假定槽内液体处于完全混合状态，由热力学可知，槽内液体温度x(t)的变化率与两种液体的温差（u(t)-x(t)）成正比，为简便计，令比例系数为1，于是有
dx(t)dtu(t)x(t) 图9.3搅拌加热器示意图 （９.1-1）
初始条件
x(0)=0 （９.1-2） 终值条件
x(1)=40 （９.1-3）
1
在1小时内散失掉的热量可用下式表示
J=[qx2(t)+ru2(t)]dt （９.1-4）
01其中q和r都是正的常数。
该问题的任务是，寻求加热液体温度的最优变化规律u*(t)，使槽内液体由初态9.1-2转移到终态9.1-3的过程中散失热量9.1-4为最小。 9.1.2最优控制问题的一般提法
根据上例可概括出，最优控制问题用数学语言来描述时应该包含以下几方面内容：
1．系统状态方程
.x=f[x(t),u(t),t] (9.1-5)
其中，x为n维状态向量，u为r维控制向量，f(.)是n维函数向量。
2．控制变量的约束条件
大多数实际控制系统中的控制变量的取值范围是受限制的，如发动机的推力，电动机的转矩等都不能超出某一极限，即 |u|≤K
在数学上，表示容许控制域为控制空间中的一个集合 uU
3．初始条件和终值条件
在最优控制中，初始条件通常是已知的，即X(t0)=X0,而终值条件则要复杂些，它可以是状态空间中一个确定的点，或状态空间中某一个点集（目标集）中的任何一个点。表示为： X(tf)=Xf
如例9.1中，tf规定为1小时。有时终值状态也可能是状态空间中一个运动的点，如用导弹拦截一个运动的目标，终值是被拦截物运动轨迹上可能的一个点，此时终值受运动轨迹约束，可表示为。
Ci(Xf、tf)=0，C=1，2，、、、，m, m 2 种不同的控制函数u(t)来实现。而要从可行的控制函数中找出一个所谓最优的控制律，首先需要的是一个评价控制效果好坏或品德优劣的性能指标。如例9.1中的热槽加温问题是把热量散失最小做为性能指标。显然，性能指标是根据控制系统本身的要求和设计者的着眼点来决定的。性能指标的一般数学形式可表示为 J=Φ[x(tf),tf]+∫ tfto L[x(t),u(t),t]dt (9.1-6) 式中第一项为终值性能指标，它表示在控制过程终了时，对系统状态变量终值的某些要求或者限制。如实际物理系统的最小稳态误差、最小定位误差等。式中的第二项为积分项性能指标，它反映在控制过程中对状态变量和控制量的某些要求和限制，在实际物理系统中，它可以表示为过程中要求各变量的综合过渡过程好、控制能量消耗最小等，如例9.1中的散热最小。 性能指标的确定，首先自然是要充分考虑物理系统本身的问题，其次也应考虑数学上的容易处理和工程上的易于实现。因此，这并不是一个纯数学问题，不存在理论指导下选定性能指标的一般方法，应根据具体情况来确定性能指标。 最优控制问题的求解，最初广泛应用变分法，但古典变分有很大局限性，它适用于状态变量及控制变量不受限制的情况，因此在工程实践中的应用受到很大限制。其后，庞特里亚金及贝尔曼分别提出了极大值原理及动态规划，为解决有限制约束的变分问题提供了有力的工具，推动了最优控制理论的发展，而如果所研究的对象为一线性系统。工程上更多采用的是易于实现和分析的二次型性能指标的线性最优控制。 9.2线性二次型最优控制 如果系统是线性的，性能指标是状态变量和控制变量的二次型函数，这样的最优控制称之为线性二次型最优控制问题。在实际工程问题中，这类问题具有特别重要的意义。其一，它可代表大量工程实际问题中提出的性能指标要求；其二，它在数学上处理比较简单，可求得最优控制的统一解析表达式。特别可贵的是，线性二次型最优控制可得到状态线性反馈的最优控制规律，易于构成闭环最优控制，从而使得其在最优控制的工程实现上具有重要的意义。 9.2.1l连续系统的二次型最优控制 1．基本概念 3 设线性定常系统为 x(t)=Ax(t)+Bu(t) y(t)=cx(t) (9.2-1) x(o)=x0 且，（C，A）是可观测的，（A，B）是可控的。 该系统的二次型最优控制问题是，寻找一个最优控制信号u(t)使得评价函数 J=..[x(t)Qx(t)+u(t)Ru(t)]dt (9.2-2) TT 0 取极小值 式中，Q为半正定对称矩阵，R为正定矩阵 式（9.2-2）所示评价函数具有两个特点： 1.该评价函数函数是状态变量和控制变量的二次型函数。 2.该评价函数是式（9.1-6）的一个特例。 此外，式（9.2-2）所示评价函数亦具有十分明确的物理含义,式中的第一项, XTQX,表示的是X到原点的一种距离,该项积分的最小,意味着要求系统状态 0从初始条件开始沿尽可能短的途径尽量向原点靠近,既考虑了过程中的偏差又考虑了终值精度。式中的第二项， 0UTRU，则表示了对控制变量u的限制,要求 控制能量代价尽可能小(因R>0,因此U→0)。 式（9.2-2）中的Q和R被称之为加权矩阵（Weight matrix），它们起着对式中两个积分项重视程度的调节钮（tuning knob）作用。即，若想重视状态的控制（如高的响应），增加Q的各元，若着重考虑控制能量消耗最小，则取较大的R值。 Q和R的选取尚无一般规律可循，一般取决于设计者的经验。常用的是所谓试行错误法（trial and error method）,即选取不同的Q、R代入计算并通过比较结果而确定。为计算方便起见，Q、R通常取为对角矩阵。 一、问题求解 对这个问题的求解可以采用多种方法，这里采用基于Liapunov 第二方法的求解方法。 4 设最优控制律为： u(t)=-kx(t) (9.2-3) 代入（9.2-1）得 x=Ax-BKx=(A-BK)x 在下面的推导中，设矩阵A-BK是稳定的，即A-BK的特征值具有负实部。 将式（9.2-3）代入式（9.2-2），得 . J=（xTQx+xTKTRKx）dt 0 =xT(Q+KTRK)xdt 0设对任意x都有 x(Q+KRK)x=T T ddt (xPx) T 式中P为正定实对称矩阵，则可导出 x(Q+KRK)x=-xPx-xPx =-x[(A-BK)P+P(A-BK)]x 由Liapunov第二方法可知，对于给定的正定矩阵Q+KTRK，如果矩阵A-BK是稳定的，则存在正定矩阵P，使得 (A-BK）TP+P（A-BK）=-（Q+KTRK） (9.2-4) 由此可求出Ｊ为： T T T T .TT .Ｊ＝ xT(Q+KTRK)xdt=-xTPx |0 0 =- xT ()p x()+xT (0)px（0） 由于假定A-BK的所有特征值都具有负实部，因此有x() →0。这样就得到 J=xT(0)Px(0) (9.2-5) 即J可由初始条件x(0)和P得到。 由于R为正定实对称矩阵，因此可将R写成 R=TT （9.2-4） T 5 式中T为一非奇异矩阵。这样式（9.2-4）可写成 （AT-KTBT）P+P（A-BK）+Q+KTTTTK=0 上式经整理还可写成 APA+[TK-（T）BP] [TK-（T）BP]-PBRBP+Q=0 J对K取极小值就要求下式对K取极小值。 xT[TK-(T)-1BTP] T[TK-(T)-1BTP]x 由于上式是非负的，其最小值为零，即 TK=（TT）-1BTP 因此可求得最优反馈增益矩阵为 K=T（T）BP=RBP （9.2-6） 最优控制u(t)为 u(t)=-Kx(t)=-R-1BTPx(t) 式（9.2-6）中的矩阵P必须满足（9.2-4）或者满足下列方程： ATP+PA-PBR-1BTP+Q=0 （9.2-7） 式（9.2-7）称为Riccati(黎卡提)方程。 二、MATLAB实现 MATLAB中提供的Lqr和lqry命令可以直接求解二次型最优控制问题（亦称最优调节器，Optimal regulator）其格式为： [K,P,E]=Lqr(A,B,Q,R) [K,P,E]=Lqry(A,B,C,D,Q,R) 式中，K为最优反馈增益矩阵；P为对应Riccati方程的正定解，E为A-BK的特征值，Lqry命令用于求解二次型调节器问题的特例，即目标函数中用输出Y来代替状态X，也就是说，目标函数为。 -1 -1 T -1 T T -1 T T -1 T -1 T J= (yTQy+uTRu)dt 0例9．2给定系统 x=Ax+Bu Y=Cx 式中 6 .0 A00102001， B0， C13100 以高的响应性能为目标，求最优控制律U=-Kx使评价指标 J=（xQx+uRu）dt 0TT 取最小化 解、取 100Q0001000; R=[0.01] 1并代入程序chp9_1.m % Program Chp9_1 % ----- Design of quadratic optimal control system ----- % ----- Design of quadratic optimal control system ----- % state matrics A,B,C,D % state matrics A,B,C,D A=[0 1 0; 0 0 1; 0 -2 -3]; B=[0; 0; 1]; C=[1 0 0]; D=0; % Weight matrices Q and R Q=[100 0 0; 0 1 0; 0 0 1]; R=[0.01]; % Obtain the optimal state feedback gain matrix K K=lqr(A,B,Q,R) % Optimal contral sysytem state matrix AA AA=A-B*K; % Unit-step response t=0:0.01:6; % Original system figure(1) [y,x,t]=step(A,B,C,D,1,t); 7 plot(t,y) grid % Optimal contral system pause figure(2) [yo,xo,to]=step(AA,B,C,D,1,t); plot(to,yo) grid 程序运行结果： 2.521.510.0120.010.0080.0060.0040.5001234560.00200123456 （a） 调节前 （b） 调节后 图 9.2 二次型最优控制 9.3离散系统的二次型最优控制 对于离散系统 x(k+1)=Ax(k)+Bu(K) Y(K)=cx(k)+Duc(k) (9.3-1) 如果考虑控制步数N→≦，则该系统的二次型最优控制问题的评价函数可表示为 J=12k0[xT(k)Qx(k)+uT(k)Ruc(k)] (9.3-2) 同样，可证明使（9.3-2）式取极小值的最优控制解可表示为： u(k)=-kx(k) (9.3-3) 式中 K=(R+BTPB)-1BTPA-Q=0 其中，P是Riccati方程 P-ATPA+ATPB（R+BTRB）-1BTPA-Q=0 8 的解 MATLAB中的 [K,P,E]=dlqr(A,B,Q,R) [K,P,E]=dlqr(A,B,C,D,Q,R) 等命令支持离散系统二次型最优控制的计算机辅助设计。 9.4状态观测器与最优观测器设计 状态反馈是现代控制中的一个重要方法，在本章的最优控制、前面讨论的极点配置以及以后的一些控制方法讨论中，都离不开状态反馈。这意味着、在控制过程中，需要对各个状态变量变化有很好的了解，而在实际控制系统中，这也意味着需要对控制变量进行在线检测。然而，实际系统中并非所有的状态变量都很容易被检测，有的状态变量可能并非物理量，有些状态变量可能很难或根本无法用物理方法去测量，有时即使状态变量可以被测量，但传感器的过多使用也可能导致系统性能的下降与成本的上升。因此，要使状态反馈方法能够很好地在工程中得以实现，需要解决这个问题。本节所讨论的状态观测器便是解决这一问题的方法之一。 所谓状态观测器是一个在物理上可实现的动力学系统，它在待观测系统的输入和输出的驱动下（这总是可以测量到的）产生一组逼近于待观测系统状态变量的输出。该动力学系统装量所输出的一组状态变量便可作为该观测系统的状态的估计值。正因为此，状态观测器也被称为状态估计器。 9．4．1 状态观测器原理 对于系统 x(t)=Ax(t)+Bu(t) (9 .4 -1) y(t)=Cx(t) 考虑一个动力学模型 （t）=Aω(t)+Ky(t)+Bu(t) ˆ(t)=Dω(t)+Hy(t) (9.4-2) x关联接系统（9.4-1）与(9.4-2)如图9.3所示 9 图9.3 状态观测器结构 如果 C1.存在一个矩阵M，并且 ˆM=MA-KC （9.4-3） A In=DM+HC （9.4-4） ˆ=MB （9.4-5） 和 BC2.A是一个稳定矩阵 ˆ(t) →x(t) (t→≦) 它将导致 x模型（9.4-2）被称为兰恩伯格观测器（Luenberger observer） 证明： 定义一个误差向量 e(t)= ω(t)-Mx(t) (9.4-6) (t)= (t)-Mx (t) 那么eˆω(t)+KCx(t)+Bˆu(t)-M[Ax(t)+Bu(t)] =Aˆω(t)-(MA-KC)x(t) =Aˆω(t)- AˆMx(t) (t)= A从而有eˆe(t) =Aˆt)e(o) (9.4-7) 所以e(t)=exp(A 10 另外，定义 ˆ(t)-x(t) (9.4-8) η（t）=x于是 η(t）=Dω(t)+HCx（t）-x（t） =Dω(t)-[In-HC]x（t） =Dω(t)-DMx（t） =De（t） ˆt)e(o) 所以η(t）=Dexp(A因为 A是一个稳定矩阵，它将导致 η（t）→0 (t→≦) ˆ(t) →x(t) (t→≦) 从而有x得证 如果让 M=D=In ,H=0 ˆ=A-KC, Bˆ=B 于是 Aˆ(t)= ω(t) xˆxˆu(t) ˆ(t)+Ky(t)+ Bˆ（t）=A从而 xˆu(t) (9.4-9) ˆ(t)+Ky(t)+ Bˆ(t)=(A-KC)x即 x式（9.4-9）所示系统即是常见的状态观测器（国外教科书上称为identify observer） 9.4.2状态观测器的设计 对式（9.4-9）进行分析可知，状态观测器设计的关健在于矩阵K的求取。一般可以通过两种方法来求取K。一是给出期望极点矢量求取，这种方法在许多教科书上均有介绍。 另一种是应用Kalman滤波进行最优估计求取，本节中介绍后一种方法。 实际系统在控制过程中，在得到状态反馈的同时，也会不可避免地受到量测噪声的干绕，所谓最优估计问题，就是指从带有量测噪声的检测状态中获得状态的最优估计。 考虑如下系统 11 (t)=Ax(t)+Bu(t)+Gv(t) xy(t)=cx(t)+w(t) 式中，v（t）为随机噪音干扰输入（System noise）；w(t)为由传感器带来的随机测量噪音（measurement noise）。进一步，假设v(t)及w(t) 为白噪声且互不相干，于是有以下统计特性： E[v(t)]=0; E[w(t)]=0 E[v(t)wT(t)]=0 E[v(t)vT(t)]=Rk E[w(t)wT(t)]=Qk 这里 Rk>0 Qk≥0 现在，对于这个包含有噪音v(t)和w(t)的系统，需要做的是根据数据u(t) ˆ(t)。 和y(t)，找出状态变量x(t)的最优估计值x 参照先前有关兰恩伯格观测器的讨论，可以假设： ˆ(t)=(A-LC)x(t)+Ly(t)+Bu(t) xˆ(o)=xo (9.4-10) x 式（9.4-10）意味着，观测器的设计取决于矩阵L的确定，这里，L被称为增益矩阵。 建立一个评价指标： T ˆ(t)][x(t)- xˆ(t)]} (9.4-11) J=E{[x(t)- x 于是，有关观测器的设计问题可转化为对（9.4-11）求极小值的问题，这个过程被称为对系统（9.4-10）的Kalman滤波。 该问题的解被求得为 L=PCRk T -1 （9.4-12） 式中，L为最优增益矩阵，或称为Kalman增益矩阵，P是下述 Riccati方程的正定解 PAT+AP-PCTRk-1CP+Qk=0 于是，通过应用最优增益矩阵L，已知输入u（t）和被测得的输出y(t)可构成一个观测器为 ˆ(t)+Ly(t)+Bu(t) ˆ(t)=[A-LC] x x 12 ˆ(t)Cyˆ(t) （9.4-13） xˆ(t)Ix 在MATLAB中，可以用命令Lqe来求解最优滤波器的参数，其格式为 [L，P，E]=Lqe（A，G，C，Qk，Rk）. 式中，L为Kalman增益矩阵，P为对应Riccati方程的解，E为估计器的闭环特征值。 命令estim用来构成连续系统Kalman滤波器，命令格式为 [Ae,Be,Ce,De]=estim(A,B,C,D,L) 根据式（9.4-13）可知，上式所求得的滤波器为含有测量噪音的传感器输出 ˆ，ˆ和状态X的估计值xY，而滤波器的作用是根据Y求出估计值y而当应用 Kalman 滤波器作状态观测器来使用时，对一些已可以用传感器检测的状态值已不必再去估计，系统的输入中也还可能有其它确定性的控制输入和外部给定输入，这时，可以应用更一般的形式的estim命令格式，如下所示。 [Ae,Be,Ce,De]=estim(A,B,C,D,L, sensors, known) 式中，Sensors是一个包含传感器输出序号的矢量，Known是包含外部输入矢量序号的矢量。有关这种情况我们将在下一节讨论，此处，我们先就比较简单的情况讨论状态估计器的设计。 例9.3给定系统 x(t)=Ax(t)+Bu(t) Y(t)=Cx(t) 式中 0A10180116210，B0，C=[0 18 –162] 100设计系统的状态估计器并进行仿真。 解：首先假设系统的输入噪音为零均值，方差为1 的高斯白噪声，传感器输 出的测量噪音为零均值，方差为0.01 的高斯白噪声，并利用randn命令形成随机噪音对系统进行仿真，然后根据Rk,Qk和命令Lqe计算最优增益矩阵L，利用estim命令形成最优估计器[Ae,Be,Ce,De]。再进行状态值和状态估计值的仿真 13 并绘出图形进行比较。 设计程序为 chp9-2.m, % Design of the estimer a0=[0 -18 162; 1 0 0; 0 1 -10]; b0=[1 0 0]';c0=[0 18 -162];d0=[0]; dimf=size(b0);noinp=dimf(2);dimc=size(c0); noout=dimc(1); q0=1;w=sqrt(q0)*randn(100,noinp); r0=0.01;nu=sqrt(r0)*randn(100,noout); t=[0:0.1:(100-1)*0.1]'; [y,x]=lsim(a0,b0,c0,d0,w,t); y=x*c0'+nu; [L,P,E]=lqe(a0,b0,c0,q0,r0);L,P,E [ae,be,ce,de]=estim(a0,b0,c0,d0,L); xh=lsim(ae,be,ce,de,y,t); subplot(2,2,1);plot(t,x(:,1),t,xh(:,2)) title('(1) x1 & xh1 versus t') subplot(2,2,2);plot(t,x(:,2),t,xh(:,3)) title('(2) x2 & xh2 versus t') subplot(2,2,3);plot(t,x(:,3),t,xh(:,4)) title('(3) x3 & xh3 versus t') xerr=x-xh(:,2:4); subplot(2,2,4),plot(t,xerr) title('(4) xerr versus t') 图9.4为状态与其估计值的对照曲线，图中（1）、（2）、（3）分别x1、x1；x2、 x2；x3、x3的对照，（4）为估计误差x=x-x的情况。 (1) x1 & xh1 versus t0.5 0.2 0(2) x2 & xh2 versus t0 -0.50.020510(3) x3 & xh3 versus t -0.200.2 5(4) xerr versus t1000 -0.020510-0.2 0510图9.4 状态与其估计值的对照曲线 14 因篇幅问题不能全部显示，请点此查看更多更全内容