92控制系统仿真及CAD 2-3 自编代码仿真

控制系统仿真及CAD——借助MATLAB自编代码仿真主讲人：王彪南京航空航天大学自动化学院2009.10控制系统MAD过程被控对象系统模型控制系统系统特性理论分析仿真分析设计理论知识演绎推理归纳总结实验数据建模验证SIMULINKMATLAB内建函数借助MATLAB自编代码实验分析球杆系统方框图与数学模型asbs1)(fsin275sgxusin7522gdtxdsincos101101lyylxx3131arctanxxyy31322231232arccoslllll23123131)()(yyxxl6224.665147.37ba直流电机+齿轮传动曲柄摇杆横梁小球数学仿真的实质系统模型数学领域微分方程传递函数状态方程控制领域频响数据数值方法求解等效离散系统仿真实质就是获取系统的动态输出响应解析解与数值解12yy已知：解析)1()(221tetyt01.02.03.04.05.0y00.43230.49080.49880.49980.5000有些微分方程难求解析解如：某些变系数微分方程有些微分方程没有解析解如：某些非线性微分方程y*00.42970.49010.49860.49980.5000数值积分法数值积分法：是求解微分方程初值问题的一种近似方法数值仿真：在仿真领域应用数值计算方法，以计算机为辅助工具，求取系统动态输出一阶微分方程（或方程组）差分递推方程集总参数系统高阶微分方程编程近似解数值积分建模1ntnt),(ytfy1),()()(1nnttnndyftyty数值积分法基本原理已知一阶微分方程00)(,),(ytyytfy)()(),(11nntttytydyfnn0tnnnQyy1积分步长：为[tn,tn+1]内f曲线下方的面积可以是向量tnntth1欧拉法与改进的欧拉法1),(nnttdyf),(nnnythfQ),(1nnnnythfyy1ntnt),(ytfyt0t1),(nnttdyf)],(),([211nnnnnytfytfhQ),(ytfyt0t)],(),([2111nnnnnnytfytfhyy欧拉法改进的欧拉法nt1nt数值积分法基本概念预估—校正法单步法计算yn+1时只需用到yn的值多步法计算yn+1时需用到ynyn-1yn-2…yn-k的值显式公式计算yn+1时所用数据均已知隐式公式(不能自启动)计算yn+1时需要用到待求量yn+1),(1nnnnythfyy)],(),([1121nnnnhnnytfytfyy),(1nnnPnythfyy)],(),([1121PnnnnhnCnytfytfyy①预估②校正nnnyhyy1hyyhyhynnnn122泰勒级数ny)(!)(!2)()()()(21nrrnnnntyrhtyhtyhtyty数值积分法原理分析（1）)(321!!3!2rnrnnnnnyrhyhyhyhyyr阶公式r=??导数??][121nnhnnyyyy数值积分法原理分析（2）),(111nnnnytfkhkyy),(),(][11212121nnnnhnnytfkytfkkkyy),(1hkyhtfnnnnnyhyy1nnnnyhyhyy!221可以利用函数f(t,y)在tn与tn+1两个离散点上的（预估）值k1、k2的线性组合，近似代替y(t)的泰勒级数中各阶导数项之和1ntnt1k2kt0tf数值积分法原理分析（3）),(),(][1212122111hkyhtfkytfkkkhyynnnnnn阶数增高可不可以呢？1ntnt1k1nt2k101t0f如何确定公式系数呢？将k2在(tn,yn)处展开成泰勒级数，并代入到yn+1的表达式中，再与y(t)在tn处的泰勒展开式的各项相比，即可确定各系数龙格—库塔法在区间[tn,tn+1]内预估多个点的函数值f(t,y)，然后用其线性组合近似代替y(t)的泰勒级数中各阶导数项之和，进而导出龙格—库塔公式1阶龙格—库塔公式即矩形公式、2阶龙格—库塔公式即梯形公式4阶龙格—库塔公式已具有较高精度，足以应付大量工程问题4阶经典龙格—库塔公式),()2,2()2,2(),()22(6342312143211hkyhtfkkhyhtfkkhyhtfkytfkkkkkhyynnnnnnnnnn4阶RK法应用举例—1阶ODE已知一阶系统G(s)，用4阶龙格—库塔法求解asbsUs)()(buayy)(),(tbuayytfy)(),(1nnnntbuayytfk)()(),(2121222hnhnhnhntbukyakytfk4阶RK法应用举例—1阶ODE（续）已知一阶系统G(s)，用4阶龙格—库塔法求解asbsUs)()(buayy)()(),(2222223hnhnhnhntbukyakytfk)()(),(334htbuhkyahkyhtfknnnn)22()()1(43216kkkknynyh4阶RK法仿真编程——1阶ODE初始化递推循环迭代绘图a、bt(1)、y(1)、u(1)h、Tnforn=1:Tn/h-1k1=f(t(n),y(n),u(n))u1=uf(t(n)+h/2)k2=f(t(n)+h/2,y(n)+h/2*k1,u1)k3=f(t(n)+h/2,y(n)+h/2*k2,u1)t(n+1)=t(n)+hu(n+1)=uf(t(n+1))k4=f(t(n)+h,y(n)+h*k3,u(n+1))y(n+1)=y(n)+h/6*(k1+2*k2+2*k3+k4)endplot(t,u,t,y)模型参数初始条件仿真参数4阶RK法应用举例—1阶ODE组已知二阶系统G(s)，用4阶龙格—库塔法求解令则)()()(assbsUsbuyayyxbuaxxbuaxxxy或12221xybuaxxxxCxBAxxyu4阶RK法应用举例—1阶ODE组（续1）已知二阶系统G(s)，用4阶龙格—库塔法求解)()()(assbsUsnnnnnxuxytfk),,,(111nnnnnnbuaxuxytfk),,,(221buaxxxy2122121122112),,,(12kxukxkytfkhnnhnhnhn1212)(),,,(2122121122222nhnnhnhnhnbukxaukxkytfk4阶RK法应用举例—1阶ODE组（续2）2222221222113),,,(12kxukxkytfkhnnhnhnhn1212)(),,,(2222221222223nhnnhnhnhnbukxaukxkytfk232313114),,,(12hkxuhkxhkyhtfknnnnn12312313224)(),,,(nnnnnnbuhkxauhkxhkyhtfk)22(6141312111kkkkhyynn)22(6242322211kkkkhxxnn如何编程？4阶RK法仿真编程—状态空间表达式已知系统状态空间表达式DuCxyBuAxx如何求解？如何处理？A、B、C、Dt、x、uh、Tnforn=1:Tn/h-1y(n)=C*x(:,n)+D*u(n)k1=A*x(:,n)+B*u(n)k2=A*(x(:,n)+h/2*k1)+B*u(t(n)+h/2))k3=A*(x(:,n)+h/2*k2)+B*u(t(n)+h/2))k4=A*(x(:,n)+h*k3)+B*u(n+1)x(:,n+1)=x(:,n)+h/6*(k1+2*k2+2*k3+k4)t(n+1)=?u(n+1)=?endplot(t,u,t,x,t,y)矩阵向量放在何处更合理？数值积分法的稳定性分析数值稳定性：用数值积分法求解微分方程，由于误差影响可能导致仿真出现不稳定的结果检验方程：稳定条件：数值稳定性与计算步长、系统特性、以及数值积分的阶数等有关yyyyrr)(nnnnnnyhyhyhyhyy4433221!4!3!21!4)(!3)(!2)(1432hhhh数值积分法的精度分析截断误差：因截去泰勒级数中高阶无穷小项而引入的误差称为r阶精度舍入误差：计算机字长有限，浮点数不能完全精确表示，因而产生误差，与字长、步长、算法等多种因素有关误差与计算步长经验步长公式)(0)(111rnnhyty截断误差舍入误差h总误差h*1min])20~5[()05.0~2.0(chTh数值积分法的速度分析主要取决于每步积分计算导数的次数矩形法1次，梯形法2次，四阶龙格—库塔法4次阶数越高，计算量越大，速度越慢在满足稳定性和精度要求的情况下，积分步长尽可能大些编程技巧也会在很大程度上影响算法的速度变步长龙格—库塔法起因：步长大了精度要降低，步长小了增加不必要的计算量，因此计算过程中合理自动调整步长非常重要基本思想：提出一个局部误差估计公式，根据估计误差取步长h，以便使估计误差始终在允许的误差范围内，从而既提高速度又提高精度具体做法：采用阶次相差一阶的两种龙格—库塔公式计算yn+1，它们的差值作为估计误差，如果大于昀大允许误差，则减小步长；如果小于昀小允许误差，则增大步长；否则不变四五阶变步长龙格—库塔法))43(2,())3(8,2())(6,3()3,3(),()4(643153142131215411kkkhyhtfkkkhyhtfkkkhyhtfkkhyhtfkytfkkkkhyynnnnnnnnnnnnRunge|Kutta|Merson)892(3054311kkkkhEn其它类型的数值积分法多步法Adams法：适用于求解非刚性微分方程；有显式公式，也有隐式公式；难以采用变步长计算面向刚性方程的算法Gear法：适用于求解刚性微分方程；有单步公式，也有多步公式；有显式公式也有隐公式；亦有变步长算法各类型方法比较：稳定性：隐式公式显式公式、多步法单步法；精度：隐式公式显式公式、多步法单步法；速度：隐式公式显式公式、多步法单步法；自启动性：隐式公式、多步法不能自启动SIMULINK中的数值积分法MATLAB中的数值积分法Runge-Kutta：ode45：4阶精度、5阶误差的变步长龙格—库塔—默森公式，单步法；大多数问题首选此法ode23：与ode45类似；但速度快些、精度低些Adams：ode113，变阶亚当姆斯—帕斯默森—梅林公式，多步法，采用预估—校正形式PECE计算；与ode45相比，容易获得高精度、计算效率高；变步长复杂Gear：ode15s，变阶吉尔法，多步法；适于求解刚性方程ode23s，与ode15s类似，但精度高些，速度慢些Ode45()应用举例functiondy=dfun(t,y)a=37.5147;b=66.6224;dy=-a*y+b*1;ode45(@dfun,[00.2],0)functiondy=dfun(t,y