1、实验三 ：加速 度指令 跟踪控 制系统 控制器 的设计 一、 实验内容 加速度指令跟踪控制系统框图如下图所示， 其中 I 为内回 路的控制器，C 为外回路 控制 器 。 、 、 分别为为 舵系 统、陀 螺仪 和加速 度计 ，设计 系统 时可选为 二 阶环节或 者 都 为 1 即可。 1 、 2 分别为从升降舵舵偏角到弹体俯仰角速率和法向加速度的传递函数，具 体环节由所给实验数据求出。 仿真时把 、 、 都当做 1 ， 1 、 2 由公式求出后， 设计控制器 C 和I，以 求系统可 以跟踪上单位阶跃信号，满足性能指标，并观察曲线得出结论。 二、 实验公式 及 数据 1 ( ) 、 2 ( ) 的
2、公式为： 1 ( ) = 3 + 2 5 3 4 2 +( 1 + 4 ) + 1 4 + 2 2 ( ) = 5 2 + 1 5 + 2 5 3 4 2 +( 1 + 4 ) + 1 4 + 2公式中各个系数所对应的数据见下表： H/km 1.067 4.526 8.21 14.288 V/m/s 546.9 609.2 701.5 880.3 a1 1.488 0.132 0.7748 0.3528 a2 104.7 91.97 76.51 46.44 a3 66.54 54.93 41.52 21.59 a4 1.296 1.126 0.9 0.514 a5 0.129 0.106 0.
3、076 0.036 三、 实验指标 实验所要求的频域和时域指标如下: 幅值裕度10dB 相角裕度45 超调量20% 上升时间0.2s 四、 设计过程 最初， 本想用 sisotool 工具进行设计， 可是试验了多次之后都失败了， 最后还是决定 用 比较熟悉的 Simiulink 框图的方法进行设计以及调试。 首先应先求出 1 ( ) 、 2 ( ) 的具体传递函数。以第一列实验数据为例。 由第一列实验数据求出 1 、 2 的具体形式为： 1 ( ) = 66.54 72.72954 2 +2.784 +106.628448C(s) G 1 (s) G d (s) G a (s) G 2 (s)
4、 n yc G g (s) I n a g zc z n y 2 ( ) = 55 0.129 2 +0.191952 72.72954 2 +2.784 +106.628448初步设想设计外环控制器 C 为 PID 控制器，内环仅采用比例控制，系统的 Simiulink 框 图如下图所示： 此时，得到的系统时域及频域图像如下： 可以看出， 在未加任何校正环节的情况下， 系统的时域响应是发散的， 很不稳定。 由开 环 Bode 图可 以看出，此时系统根本就没有穿越 0db 线，没有剪切频率，高频也没有衰减， 导致系统还没有进行工作就已经不稳定了。 因此系统需要添加控制器， 通过控制器的校正使系
5、统能在工作频率内处于稳定状态。 初 步设计加入 PID 控制器。 加 PID 控制器的原因有： （1 ） 在 Simiulink 中，PID 控制器调试非常的方便； （2 ） P 控制器可以使系统的开环 Bode 图向上平移，增加系统的带宽； （3 ） I 控制器可以提高系统的型别，起到消除稳态误差的作用； （4 ） D 控制器起 到了增加阻尼的作用，可以提高系统的相对稳定性。 因此，PID 控制器虽然结构简单，但是作用不容小觑。 经过反复调试 PID 的参数，初步定为： P = 0.06 I = 1.1 D = 0.01 在此参数下得到系统的的单位阶跃响应图像如下图， 可见此时系统的响应速度
6、很快， 但 是超调量很大， 振荡次数很多， 向下的过冲也很大， 明显不符合实验指标要求。 根据现在 的 结果，需要增加环节滤去系统中比较高频的分量，让系统的变化更加平稳光滑。 因此， 为了滤去系统中比较高频的分量经过反复尝试， 最后在内环的前向通道加了一个 惯性环节： 加入这个惯性环节， 相当于让系统的带宽降低， 使高频段被及时的衰减， 此时， 系 统的 时域响应和频域响应都有很大的好转。 最终，系统仿真的 Simiulink 框图如下 图所示： 图中 这个地方原来应为 ，在 55 左 右， 变为 20 就相 当于加入了一个比例环节 20 55 。 其中，内环的参数及结构为： PID 的具体参
7、数为： 五、 试验结果 及分 析 试验仿真得 出的时域跟 踪结果有过 载输出，与 事先预想的 一样，具有 非最小相位 特性 。 产生此现象的物理原因为： 若使导弹产生正方向的法向过载， 需要舵偏角负偏， 产生正的攻 角， 使导弹抬头。 而舵偏转的速度非常快， 远快于导弹攻角变化的速度。 也就是舵已经很 快 地负偏， 但是攻角还没有建立起来。 因为舵负偏， 产生了一个负的法向过载。 所以阶跃响 应 是先向反方向变化。 当攻角逐渐地建立起来， 攻角产生的正过载， 抵消了舵负偏产生的负过 载。攻角逐渐地变大，最终到达所需位置，跟踪上过载指令。 由实验结果可以看出，此时系统的性能指标为： 超调量在 2
8、0% 左右 调整时间不到 2s 上升时间为 0.5s 实验结果比 实验要求的 指标要大一 些，但基本 符合要求。 时域响应有 明显的过载 现象 。 与在前没有加入惯性环节相比， 系统的性能有了很大的提升。 时域指标基本满足之后， 观察 频域指标。 系统的频域响应如下， 下图分别为开环系统的幅频及相频特性以及闭环系统的幅频及相 频特性。由开环系统的幅频及相频特性可以看出系统的： 剪切频率为 3.09rad/s 幅值裕度为 8.52dB 相角裕度为 48 基本符合系统要求的频域性能指标，因此设计基本成功。 六、 实验结论 本次加速度指令控制系统控制器设计实验，内环控制器选用了一个一阶惯性环节 ，
9、外环控制器 C 选用了 PID 控制器， 其中P = 0.06 、I = 1.1 、D = 0.01 ， 在两者 的作用下， 系统基本可以达到事先所要求的指标。 之所以在内环加入了一个一阶惯性环节的 原因其实我也说不太清楚， 本希望内环只用比例环节就可以实现， 可是在怎么调试都不行之 后， 就觉得应该在内环加点什么， 为了衰减高频分量， 就加入了一阶惯性。 其实， 这样调试 的步骤是有待商榷的， 因为系统调试的一般步骤为先调试内环， 再调试外环， 所以下次这种 任性的行为还是应该避免。 本次实验， 虽然没有学会 sisotool 的使用， 但是让我在 Simiulink 中的 PID 调试更加
10、的熟 练， 对比例、 积分、 微分环节各自对系统的作用更加了解， 并积累了调试系统的经验， 可以 作为以后学习和工作的宝贵财富。 实验四 ： 弹道 仿真 一、 实验内容 目标的运动情况自己假定， 选择合适的比例导引系数， 利用四阶龙格库塔求解出仿真结 果，绘出导弹与目标的运动轨迹。 二、 实验公式 及 数据 导弹的动力学和运动学方程同实验 2，如 式 (2)(3)(4)(5)(6)(7)所示 ， 目标的动力学方程为： 0 = t v (1) cos = t t yt t v ag (2) cos = t t vt tz va (3) 目标的运动学方程为： cos cos = t t t vt
11、xv (4) sin = tt t yv (5) cos sin = t t t vt zv (6) 比例导引律： = y a KVq (7) = z a KVq (8) 其中， 22 arctan = + r rr y q xz(9) arctan = r r z q x(10) 目标相对导弹的运动方程： = rt x xx (11) = rt y yy (12) = rt z zz (13) 其中， = rt x xx (14) = rt y yy (15) = rt z zz (16) 三、 实验内容 1. Simulink 中的结构图 2. 各部分 s 函数的编制 2.1 导弹动力 学
12、部分（MissileDynamic ） function sys,x0,str,ts,simStateCompliance = MissileDynamic(t,x,u,flag) switch flag, case 0, sys,x0,str,ts,simStateCompliance=mdlInitializeSizes; case 1, sys=mdlDerivatives(t,x,u); case 2, sys=mdlUpdate(t,x,u); case 3, sys=mdlOutputs(t,x,u); case 4, sys=mdlGetTimeOfNextVarHit(t,x,
13、u); case 9, sys=mdlTerminate(t,x,u); otherwise DAStudio.error(Simulink:blocks:unhandledFlag, num2str(flag); end function sys,x0,str,ts,simStateCompliance=mdlInitializeSizes sizes = simsizes; sizes.NumContStates = 3; sizes.NumDiscStates = 0; sizes.NumOutputs = 3; sizes.NumInputs = 2; sizes.DirFeedthr
14、ough = 1; sizes.NumSampleTimes = 1; sys = simsizes(sizes); x0 = -40/180*pi 0 450; str = ; ts = 0 0; simStateCompliance = UnknownSimState; function sys=mdlDerivatives(t,x,u) v=x(3); g=9.8; sys(1)=(u(1)-cos(x(1)*g/v; %sita sys(2)=-u(2)*g/(v*cos(x(1); %fc sys(3)=-g*sin(x(1); %v function sys=mdlUpdate(t
15、,x,u) sys = ; function sys=mdlOutputs(t,x,u) sys(1)=x(1); %sita sys(2)=x(2); %fc sys(3)=x(3); %v 2.2 导弹运动 学部分（MissileMovtion ） function sys,x0,str,ts,simStateCompliance = MissileMovtion(t,x,u,flag) switch flag, case 0, sys,x0,str,ts,simStateCompliance=mdlInitializeSizes; case 1, sys=mdlDerivatives(t
16、,x,u); case 2, sys=mdlUpdate(t,x,u); case 3, sys=mdlOutputs(t,x,u); case 4, sys=mdlGetTimeOfNextVarHit(t,x,u); case 9, sys=mdlTerminate(t,x,u); otherwise DAStudio.error(Simulink:blocks:unhandledFlag, num2str(flag); end function sys,x0,str,ts,simStateCompliance=mdlInitializeSizes sizes = simsizes; si
17、zes.NumContStates = 0; sizes.NumDiscStates = 0; sizes.NumOutputs = 3; sizes.NumInputs = 3; sizes.DirFeedthrough = 1; sizes.NumSampleTimes = 1; sys = simsizes(sizes); x0 = ; str = ; ts = 0 0; simStateCompliance = UnknownSimState; function sys=mdlDerivatives(t,x,u) sys = ; function sys=mdlUpdate(t,x,u
18、) sys = ; function sys=mdlOutputs(t,x,u) v=u(3); sys(1)=v*cos(u(1)*cos(u(2); sys(2)=v*sin(u(1); sys(3)=-v*cos(u(1)*sin(u(2); 2.3 目标运动 模型（TargetMovtion ） function sys,x0,str,ts,simStateCompliance = TargetMovtion(t,x,u,flag) switch flag, case 0, sys,x0,str,ts,simStateCompliance=mdlInitializeSizes; cas
19、e 1, sys=mdlDerivatives(t,x,u); case 2, sys=mdlUpdate(t,x,u); case 3, sys=mdlOutputs(t,x,u); case 4, sys=mdlGetTimeOfNextVarHit(t,x,u); case 9, sys=mdlTerminate(t,x,u); otherwise DAStudio.error(Simulink:blocks:unhandledFlag, num2str(flag); end function sys,x0,str,ts,simStateCompliance=mdlInitializeS
20、izes sizes = simsizes; sizes.NumContStates = 3; sizes.NumDiscStates = 0; sizes.NumOutputs = 3; sizes.NumInputs = 2; sizes.DirFeedthrough = 1; sizes.NumSampleTimes = 1; sys = simsizes(sizes); x0 = 0 0 272; str = ; ts = 0 0; simStateCompliance = UnknownSimState; function sys=mdlDerivatives(t,x,u) g=9.
21、8; nyt=u(1); nzt=u(2); sys(1)=g*(nyt-cos(x(1)/x(3); %sita sys(2)=-g*nzt/(x(3)*cos(x(1); %fc sys(3)=-g*sin(x(1); %v function sys=mdlUpdate(t,x,u) sys = ; function sys=mdlOutputs(t,x,u) sys(1)=x(3)*cos(x(1)*cos(x(2); sys(2)=x(3)*sin(x(1); sys(3)=-x(3)*cos(x(1)*sin(x(2); 2.4 导引头测量（Measuring ） function
22、sys,x0,str,ts,simStateCompliance = Measuring(t,x,u,flag) switch flag, case 0, sys,x0,str,ts,simStateCompliance=mdlInitializeSizes; case 1, sys=mdlDerivatives(t,x,u); case 2, sys=mdlUpdate(t,x,u); case 3, sys=mdlOutputs(t,x,u); case 4, sys=mdlGetTimeOfNextVarHit(t,x,u); case 9, sys=mdlTerminate(t,x,u
23、); otherwise DAStudio.error(Simulink:blocks:unhandledFlag, num2str(flag); end function sys,x0,str,ts,simStateCompliance=mdlInitializeSizes sizes = simsizes; sizes.NumContStates = 0; sizes.NumDiscStates = 0; sizes.NumOutputs = 3; sizes.NumInputs = 8; sizes.DirFeedthrough = 1; sizes.NumSampleTimes = 1
24、; sys = simsizes(sizes); x0 = ; str = ; ts = 0 0; simStateCompliance = UnknownSimState; function sys=mdlDerivatives(t,x,u) sys = ; function sys=mdlUpdate(t,x,u) sys = ; function sys=mdlOutputs(t,x,u) %r=zeros(3,1); Vr=u(1);u(2);u(3); r=u(4);u(5);u(6); O=cross(r,Vr)/(u(4)2+(u(5)2+(u(6)2); s=u(7); f=u
25、(8); A=cos(s)*cos(f) sin(s) -cos(s)*sin(f); -sin(s)*cos(f) cos(s) sin(s)*sin(f); sin(f) 0 cos(f); de=A*O; dr=dot(r,Vr)/sqrt(u(4)2+(u(5)2+(u(6)2); sys(1)=dr; %r sys(2)=de(2); %deew sys(3)=de(3); %deez 2.5 测量误差（MeasureError） function sys,x0,str,ts,simStateCompliance = MeasureError(t,x,u,flag) switch f
26、lag, case 0, sys,x0,str,ts,simStateCompliance=mdlInitializeSizes; case 1, sys=mdlDerivatives(t,x,u); case 2, sys=mdlUpdate(t,x,u); case 3, sys=mdlOutputs(t,x,u); case 4, sys=mdlGetTimeOfNextVarHit(t,x,u); case 9, sys=mdlTerminate(t,x,u); otherwise DAStudio.error(Simulink:blocks:unhandledFlag, num2st
27、r(flag); end function sys,x0,str,ts,simStateCompliance=mdlInitializeSizes sizes = simsizes; sizes.NumContStates = 0; sizes.NumDiscStates = 0; sizes.NumOutputs = 3; sizes.NumInputs = 3; sizes.DirFeedthrough = 1; sizes.NumSampleTimes = 1; sys = simsizes(sizes); x0 = ; str = ; ts = 0 0; simStateComplia
28、nce = UnknownSimState; function sys=mdlDerivatives(t,x,u) sys = ; function sys=mdlUpdate(t,x,u) sys = ; function sys=mdlOutputs(t,x,u) sys(1)=u(1)+0.001; sys(2)=u(2)+0.002; sys(3)=u(3)+0.0005; % 2.6 导引律（GuideRule ） function sys,x0,str,ts,simStateCompliance = GuideRule(t,x,u,flag) switch flag, case 0
29、, sys,x0,str,ts,simStateCompliance=mdlInitializeSizes; case 1, sys=mdlDerivatives(t,x,u); case 2, sys=mdlUpdate(t,x,u); case 3, sys=mdlOutputs(t,x,u); case 4, sys=mdlGetTimeOfNextVarHit(t,x,u); case 9, sys=mdlTerminate(t,x,u); otherwise DAStudio.error(Simulink:blocks:unhandledFlag, num2str(flag); en
30、d function sys,x0,str,ts,simStateCompliance=mdlInitializeSizes sizes = simsizes; sizes.NumContStates = 2; sizes.NumDiscStates = 0; sizes.NumOutputs = 2; sizes.NumInputs = 4; sizes.DirFeedthrough = 0; sizes.NumSampleTimes = 1; sys = simsizes(sizes); x0 = 0 0; str = ; ts = 0 0; simStateCompliance = Un
31、knownSimState; function sys=mdlDerivatives(t,x,u) Tg=0.01; % g=9.8; K=3; %k=2:4 deew=u(2); deez=u(3); dr=u(1); sys(1)=(K*abs(dr)*deez/g+cos(u(4)-x(1)/Tg; sys(2)=(-K*abs(dr)*deew*cos(u(4)/g-x(2)/Tg; function sys=mdlUpdate(t,x,u) sys = ; function sys=mdlOutputs(t,x,u) sys(1)=x(1); sys(2)=x(2); 3. 仿真结果
32、及分析 对照初始参数：导弹：v m =450m/s =-40 =0 x m =0 y m =1000m z m =0 目标：v =272m/s =0 =0 x=1000m y=0 z=0 K 值设定为 3 仿真结果图如下： 导弹过载 n yc 与 n zc 变化如 下： 4. 调整参数 再次运行 K=1 时（弹道 发散，不能打到目标） ： K=2 ： K=4 ： K=6 ： K=10： 可见，随着比例系数 k 的增加导弹弹道趋于收敛，导引弹道更加平直，击中目标时需 用过载较小，但 k 过大 时初始段需用过载过大，实际 中 导 弹 将 不 能 按 照 比 例 导引弹道飞行， 并且 k 过大 影响
33、制导系统稳定性， 外界干扰将对导弹飞行影响增大， 而 k 过小时过载稳定性 由图可见不好，且弹道较弯甚至发散以致导弹不能击中目标，故 k 一 般应取 2-4 之间 。 四、 总结 经过这次仿真实践，我们对比例导引法有了进一步认识，熟悉了导弹运动模型，掌握 了 MATLAB 编程及用 Simulink 仿真的 一些基本知识。 其他的感受最深的是在编程过程中， 不 时会有一些小错误出现，不过收获很大。 飞行器制导与控制 实验报告 班 级： 1204201 班 学 号： 112011412 姓 名： 白赫 题 目：加速度指令控制系统设计与弹道仿真 指导教师： 周荻 2015 年 12 月 17 日
