用MATLAB解最优控制问题及应用实例学习教案

1、会计学1用用MATLAB解最优控制问题及应用解最优控制问题及应用(yngyng)实例实例第一页，共128页。第1页/共128页第二页，共128页。第2页/共128页第三页，共128页。第3页/共128页第四页，共128页。DuCxyBuAxx 1, 系统模型的建立 系统的状态方程为：第4页/共128页第五页，共128页。,;,;,;,;,;,;,;,;,212222111211212222111211212222111211212222111211qpqqppqnqqnnnpnnppnnnnnndddddddddDcccccccccCbbbbbbbbbBaaaaaaaaaA第5页/共128页第

2、六页，共128页。传递函数的零极点模型为：)()()()()(2121nmpspspszszszsKsG在MATLAB中可以采用如下语句将零极点模型输入到工作(gngzu)空间:;2121nmpppPzzzZKKGainzpk(Z,P,KGain)第6页/共128页第七页，共128页。传递函数模型在更一般(ybn)的情况下，可以表示为复数变量s的有理函数形式：nnnnnmmmmasasasasbsbsbsbsG122111121)(第7页/共128页第八页，共128页。在MATLAB中可以采用如下语句将以上的传递函数模型输入到工作空间:G=tf(num,den);, 1 ;,121121nnm

3、maaaadenbbbbnum第8页/共128页第九页，共128页。第9页/共128页第十页，共128页。第10页/共128页第十一页，共128页。第11页/共128页第十二页，共128页。第12页/共128页第十三页，共128页。第13页/共128页第十四页，共128页。第14页/共128页第十五页，共128页。)()()()()(tUtBtXtAtX)()()(tXtCtY( )X tn( )U tm)(tYl第15页/共128页第十六页，共128页。寻找(xnzho)最优控制，使下面的性能指标最小fttTTffTdttUtRtUtetQtetPeteuJ0)()()()()()(21)(

4、)(21)(l其中， 是 对称半正定常数阵， 是 对称半正定阵, 是 对称正定阵。llll)(tQ)(tRmmP第16页/共128页第十七页，共128页。1( )( ) ( )( ) ( )( ) ( )( )( ) ( )( )TTK tK t AtA t K tK t Bt R t B t K tQt第17页/共128页第十八页，共128页。可以看出，上述的黎卡提矩阵微分方程求解起来非常困难，所以我们往往求出其稳态解。例如目标函数中指定终止时间可以设置成 ，这样可以保证系统状态渐进的趋近于零值，这样可以得出矩阵趋近于常值矩阵 ，且 ，这样上述黎卡提矩阵微分方程可以简化成为：ft)(tK0)

5、(tK1( ) ( )( )( )( ) ( )( )( )( )( )0TTK t A tAt K tK t B t Rt Bt K tQ t第18页/共128页第十九页，共128页。这个方程称为代数黎卡提方程。代数黎卡提方程的求解非常简单，并且其求解只涉及到矩阵(j zhn)运算，所以非常适合使用MATLAB来求解。第19页/共128页第二十页，共128页。求解代数黎卡提方程的算法有很多，下面我们介绍一种简单的迭代算法来解该方程，令 ，则可以写出下面的迭代公式00 11TTTTiiiiiEEEGWGGHEGWQ第20页/共128页第二十一页，共128页。AIAIE 1BAIG12BAIQA

6、IBRHTT11)(BAIQW1第21页/共128页第二十二页，共128页。如果 收敛于一个常数矩阵，即 ，则可以得出代数黎卡提方程的解为：上面的迭代算法可以用MATLAB来实现：1iii 11112AIAIPiT第22页/共128页第二十三页，共128页。第23页/共128页第二十四页，共128页。第24页/共128页第二十五页，共128页。第25页/共128页第二十六页，共128页。这里的求解是建立在MATLAB的控制系统工具箱中给出的一个基于Schur变换的黎卡提方程求解函数are()基础上的，该函数的调用(dioyng)格式为： X=are(M,T,V)第26页/共128页第二十七页，

7、共128页。其中， 矩阵满足下列代数黎卡提方程，are是Algebraic Riccati Equation的缩写。对比前面给出的黎卡提方程，可以容易得出VTM,0VXTXXMMXTAMTBBRT1QV第27页/共128页第二十八页，共128页。第28页/共128页第二十九页，共128页。第29页/共128页第三十页，共128页。线性系统为： ，其目标函数是：uxx10351006667. 1 10020020050021dtuuxxJTT确定(qudng)最优控制。第30页/共128页第三十一页，共128页。第31页/共128页第三十二页，共128页。第32页/共128页第三十三页，共128

8、页。第33页/共128页第三十四页，共128页。第34页/共128页第三十五页，共128页。第35页/共128页第三十六页，共128页。第36页/共128页第三十七页，共128页。)(6.7496)(13.0276)(21*txtxtu第37页/共128页第三十八页，共128页。第38页/共128页第三十九页，共128页。第39页/共128页第四十页，共128页。图图12-1 最优控制曲线与最优最优控制曲线与最优状态曲线状态曲线第40页/共128页第四十一页，共128页。第41页/共128页第四十二页，共128页。第42页/共128页第四十三页，共128页。第43页/共128页第四十四页，共1

9、28页。)(6.7775)(13.0608)(21*txtxtu第44页/共128页第四十五页，共128页。无人飞行器的最优高度控制(kngzh)，飞行器的控制(kngzh)方程如下 )(2/100)()()(2/100100010)()()(tuthththththth 第45页/共128页第四十六页，共128页。 是飞行器的高度； 是油门输入；设计控制律使得如下指标最小)(th)(tudttuththththththtutxJ02)(2)()()(000000001)(),(),(21)(),( 初始状态 。绘制系统状态与控制输入，对如下给定的 矩阵进行仿真分析.Tththth0 , 0

10、,10)(),()( ，,QR第46页/共128页第四十七页，共128页。100000 ,R2000Q2000R,000000001Q2,0000000010RQ10001000 ,2000QR第47页/共128页第四十八页，共128页。第48页/共128页第四十九页，共128页。第49页/共128页第五十页，共128页。图图12-2 12-2 状态响应曲线状态响应曲线(qxin)(qxin)及控制输入响应曲线及控制输入响应曲线(qxin)(qxin)第50页/共128页第五十一页，共128页。第51页/共128页第五十二页，共128页。第52页/共128页第五十三页，共128页。图图12-3

11、 12-3 状态响应曲线及控制状态响应曲线及控制(kngzh)(kngzh)输入响应曲线输入响应曲线第53页/共128页第五十四页，共128页。第54页/共128页第五十五页，共128页。第55页/共128页第五十六页，共128页。图图12-4 12-4 状态响应曲线及控制状态响应曲线及控制(kngzh)(kngzh)输入响应曲线输入响应曲线第56页/共128页第五十七页，共128页。第57页/共128页第五十八页，共128页。第58页/共128页第五十九页，共128页。图图12-5 12-5 状态响应曲线及控制状态响应曲线及控制(kngzh)(kngzh)输入响应曲线输入响应曲线第59页/共

12、128页第六十页，共128页。由1），2），3），4）可分析(fnx)如下：图12-3与图12-2相比，当Q不变，R增大时，各相应曲线达到(d do)稳态所需时间增长，即响应变慢；但波动幅值变小，反馈矩阵变小；图12-4与图12-2和图12-3相比，当Q对角线上第1个元素增大时，各相应曲线达到稳态所需时间(shjin)变短，即响应快；但波动幅值变大，反馈矩阵增大；第60页/共128页第六十一页，共128页。由图12-5可知，当Q对角线上第2个元素(yun s)增大时，状态x1，x2曲线达到稳态所需时间较长，即响应较慢，平缓的趋于零；状态x3，控制输入u达到稳态所需时间短，即响应快；状态x2，x

13、3波动幅值较小，比图12-2和图12-4小，比图12-3稍大，控制输入u波动幅值比图12-2和图12-4小，比图12-3大；反馈矩阵最大。第61页/共128页第六十二页，共128页。 综上所述可得结论：Q=diag（1，0，0），R=2时，系统各方面响应较好。 矩阵Q变大时，反馈矩阵变大；当Q的对角线上第1个元素(yun s)变大时，各曲线波动幅值变大，达到稳态所需时间变短； 第62页/共128页第六十三页，共128页。 当Q的对角线上第2个元素变大时，各曲线(qxin)波动幅值变小；达到稳态所需时间，状态x1，x2增长，状态x3，控制输入u变短； 当R变大时，反馈矩阵变小；各曲线(qxin)

14、波动幅值变小；达到稳态所需时间变长。 所以根据实际的系统允许，我们应该适当选择Q和R。第63页/共128页第六十四页，共128页。第64页/共128页第六十五页，共128页。第65页/共128页第六十六页，共128页。第66页/共128页第六十七页，共128页。第67页/共128页第六十八页，共128页。第68页/共128页第六十九页，共128页。第69页/共128页第七十页，共128页。第70页/共128页第七十一页，共128页。第71页/共128页第七十二页，共128页。第72页/共128页第七十三页，共128页。第73页/共128页第七十四页，共128页。第74页/共128页第七十五页，

15、共128页。第75页/共128页第七十六页，共128页。第76页/共128页第七十七页，共128页。在现有的自寻的导弹中，大都采用比例导引法。假设导弹和目标在同一平面内运动，按比例导引制导律，假设导弹的速度向量的旋转角速度 垂直于瞬时的弹目视线，并且正比于导弹与目标之间的视线角速率 ，假设目标的法向加速度为零，那么可得： (12-1)q qN第77页/共128页第七十八页，共128页。其中， 为导弹的速度与基准方向的夹角， 为导弹与目标连线与基准方向的夹角，称为视线角， 是视线角速率， 是比例常数，称为导航比，通常为36。比例导引的实质是使导弹向着 减小的方向运动，抑制视线旋转，也就是使导弹的

16、相对速度对准目标，保证导弹向着前置碰撞点飞行。qq Nq 第78页/共128页第七十九页，共128页。第79页/共128页第八十页，共128页。第80页/共128页第八十一页，共128页。假设导弹和目标在同一平面内运动，如图12-6所示。选 为固定坐标。导弹速度向量 与 轴成 角，目标速度向量 为与 轴成 角。导弹与目标的连线 与轴 成 角。假定导弹以尾追的方式攻击目标。坐标轴 和 的方向可以任意选择，使 和 都比较小。再假定导弹和目标均匀速飞行，也就是说 和 均为恒值。使用相对坐标状态变量，设 为导弹与目标在 轴方向上的距离偏差， 为导弹与目标在 轴方向上的距离偏差，即 (12-2)oxyM

17、VoyTVoyTMToyqoxoyT,qMVTVxoxyoyMTMTyyyxxx第81页/共128页第八十二页，共128页。xxyyMxTxMyTyMVTV0qTRTMH图图12-6 12-6 导弹和目标运动几何导弹和目标运动几何(j h)(j h)关系图关系图第82页/共128页第八十三页，共128页。假定 和 比较小，因此，则1cos, 1cos,sin,sinTTTT将上式对t求导，并根据导弹和目标的关系（如图12-6所示）可得 (12-3) coscossinsinMTTMTMTTMTVVyyyVVxxxMTMTMTTMTVVyyyVVxxx (12-4)第83页/共128页第八十四页

18、，共128页。以 表示 ， 表示 （即 ），则 (12-5) (12-6)1xx2xx 1x 21xx MTTVVxx2式中 表示目标的横向加速度， 表示导弹横向加速度，分别以 和 表示，那么 (12-7)TTVMVTaMaMTaax2第84页/共128页第八十五页，共128页。导弹的横向加速度 为一控制量。一般将控制信号加给舵机，舵面偏转后产生弹体攻角 ，而后产生横向加速度 。如果忽略舵机和弹体的惯性，而且假设控制量的单位与加速度单位相同，则可用控制量 来表示 ，也就是令 (12-8) 所以(12-7)式为： (12-9)MaMauMaMauuaxT2第85页/共128页第八十六页，共128

19、页。这样可得导弹运动状态方程为： (12-10) (12-11)21xx Taux2第86页/共128页第八十七页，共128页。可写成矩阵的形式： (12-12)式中， ， ， ， 。 (12-13)如果不考虑目标的机动，即 ，则在这种情况下，式(12-12)变成： (12-14)TDaBuAXX21xxX0010A10B10D0TaBuAXX第87页/共128页第八十八页，共128页。下面来考虑(12-4)式，该式可写成 (12-15) 其中 表示导弹相对目标的接近速度。由于 的值都比较小， 可近似表示导弹与目标之间的距离。设 为导弹与目标的遭遇时刻(即导弹与目标相碰撞或两者之间的距离为最短

20、的时刻)，则在某一瞬时 ，导弹与目标的距离 可近似用下式表示： (12-16)(TMVVyCTMVVVT,和qyftty)()()(ttVttVVtyfCfTM第88页/共128页第八十九页，共128页。又考虑到对于导弹制导来说，最基本的要求是脱靶量越小越好，因此，应该选择最优控制量 ，使得下面的指标函数为最小。 (12-17)u22)()()()(fMfTfMfTtytytxtxJ第89页/共128页第九十页，共128页。然而，当要求一个反馈形式的控制时，按上式列出的问题很难求解。所以我们以 时刻，即 时 的值 作为脱靶量，要求 值越小越好。另外，由于舵偏角受到限制，导弹结构能够承受的最大载

21、荷也受到限制，所以控制信号 也应该受到限制。 ftt 0)()(ffCfttVty)(1ftx)(1ftxu第90页/共128页第九十一页，共128页。因此，我们选择以下形式的二次型指标函数： (12-18)式中 ， 。 (12-19)即 (12-20)给定初始条件 ，应用最优控制理论，可以求出使 为最小的 。dtRuuQXXtCXtXJfttTTffT021212100ccC0000QdtRutCXtXJfttffT022121 0tXJu第91页/共128页第九十二页，共128页。第92页/共128页第九十三页，共128页。 (12-21)式中 满足下列黎卡提矩阵微分方程 (12-22)

22、的终端条件为 (12-23)因此求解线性二次型问题的关键是求解黎卡提矩阵微分方程。 KXBRtuT1*QKBKBRKAKAKTT1KKCtKf第93页/共128页第九十四页，共128页。当不考虑弹体惯性时，而且假定目标不机动，即，导弹运动状态方程为 (12-24)BuAXX第94页/共128页第九十五页，共128页。指标函数为 (12-25)式中 ， , ， 。 dtRuuQXXtCXtXJfttTTffT021210010A10B2100ccC0000Q第95页/共128页第九十六页，共128页。给出 时刻， 的初值 ，采用极小值原理可求得最优控制 为 (12-26)0tt 21xx 和 0

23、201txtx和 tu* 2321 21 211212*34211 22231312fffffffcc ttcc ttc ttxcc ttxRRu tc ttc ttcc ttRRRR 第96页/共128页第九十七页，共128页。在指标函数中，如不考虑导弹的相对运动速度 项，则可令 。 变成 (12-27)以 除上式的分子和分母，得 (12-28)2x02c tu* RttcRxttcxttctufff313122111*1c 31221*333ttcRxttxtttufff第97页/共128页第九十八页，共128页。为了使脱靶量为最小，应选取 ，则 (12-29)根据图12-6可得 1c t

24、txttxtuff221*3ttVxyxtgqfc11第98页/共128页第九十九页，共128页。当 比较小时， ，则 (12-30) (12-31)qqtgq ttVxqfc1ttxttxVttVxttxqffcfcf2212111第99页/共128页第一百页，共128页。将上式代入(12-29)式，可得 (12-32)在上式中， 的单位是加速度的单位 。把 与导弹速度向量 的旋转角速度 联系起来，则有 (12-33) qVtuc3*u2秒米DVqVVVuMcM3u第100页/共128页第一百零一页，共128页。从(12-32)和(12-33)式可以看出，当不考虑弹体惯性时，最优导引规律就是

25、比例导引，其导航比为 。这证明了比例导引是一种很好的导引方法。最优导引规律的形成可用图12-7来表示。McVV3第101页/共128页第一百零二页，共128页。cV3s1s1ttVfc121ttVfcq*u1x2x图图12-7 12-7 最优导引最优导引(do yn)(do yn)方框图方框图第102页/共128页第一百零三页，共128页。xxyyMxTxMyTyMVTV0HqTTMRTaMaMTLHE图图12-8 12-8 最优导引攻击几何最优导引攻击几何(j h)(j h)平面平面第103页/共128页第一百零四页，共128页。最优导引攻击几何关系如图12-8所示，在这里讨论的目标和导弹均

26、认为是二维拦截几何平面上的质点，分别以速度 和 运动。导弹的初始位置为相对坐标系的参考点，导弹初始速度矢量指向目标的初始位置， 为导弹的指令（垂直于视线）。TVMVMa第104页/共128页第一百零五页，共128页。其中： (12-34) (12-35) (12-36) 为目标速度在 轴上的分解， 是目标的角度。导弹和目标之间的接近速度为： (12-37)TTTVaTTTyVVcosTTTxVVsinTcTMVR ,TxTyVV, x y第105页/共128页第一百零六页，共128页。目标的速度分量可由其位置变化得到： (12-38)同样地，我们可以得到导弹的位置和速度的微分方程： ， (12

27、-39) ， (12-40)TxTxTyTyVRVR,MxMxaVMyMyaVMxMxVRMyMyVR第106页/共128页第一百零七页，共128页。上面几式中的下标x，y分别表示在x和y轴上的分量。 是导弹在地球坐标系的加速度分量。为了得到导弹的加速度分量，我们必须得到弹目的相对位移： (12-41) (12-42)MyMxaa,MxTxTMxRRRMyTyTMyRRR第107页/共128页第一百零八页，共128页。TMyTMxRRq1tanMxTxTMxVVVMyTyTMyVVV从图12-8中，根据三角关系我们可以得到视线角： (12-43)如果定义地球坐标系的速度分量为： (12-44) (12-45)第108页/共128页第一百零九页，共128页。我们可以根据视线角的公式求导后得到视线角速率： (12-46) (12-47)所以我们不难得出弹目的接近速度为： (12-48)TMTMy