【利用母船和牵引风箬在飞行中回收微型飞行器（MAVs）的方法】使用高斯原理推导了电缆-风箬系统的动力学模型（Matlab、Simulink仿真实现）

   💥💥💞💞欢迎来到本博客❤️❤️💥💥

🏆博主优势：🌞🌞🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。

⛳️座右铭：行百里者，半于九十。

 ⛳️赠与读者

👨‍💻做科研，涉及到一个深在的思想系统，需要科研者逻辑缜密，踏实认真，但是不能只是努力，很多时候借力比努力更重要，然后还要有仰望星空的创新点和启发点。当哲学课上老师问你什么是科学，什么是电的时候，不要觉得这些问题搞笑。哲学是科学之母，哲学就是追究终极问题，寻找那些不言自明只有小孩子会问的但是你却回答不出来的问题。建议读者按目录次序逐一浏览，免得骤然跌入幽暗的迷宫找不到来时的路，它不足为你揭示全部问题的答案，但若能让人胸中升起一朵朵疑云，也未尝不会酿成晚霞斑斓的别一番景致，万一它居然给你带来了一场精神世界的苦雨，那就借机洗刷一下原来存放在那儿的“躺平”上的尘埃吧。

     或许，雨过云收，神驰的天地更清朗.......🔎🔎🔎

💥1 概述

参考文献1：

摘要：本文提出了一种利用母船和牵引风箬在飞行中回收微型飞行器（MAVs）的方法。基于高斯原则提出了母船-电缆-风箬系统动力学建模的方法。利用系统的微分平整性质计算从期望的风箬轨道到母船轨迹的方法，并提出了基于李亚普诺夫原理的控制器，实现了准确的母船轨迹跟踪。此外，还描述了用于风箬的基于阻力的控制器。讨论了使MAV估计和跟踪风箬轨道的方法。通过仿真和飞行结果展示了建模和控制方法。

关键词： 无人空中车辆；MAVs；微型飞行器；动态建模；轨迹跟踪；多车辆；空中回收；牵引体系统；反步法；轨道估计；轨道跟踪；自主系统。

参考文献2： 

摘要-本文提出了一种利用大型母船和回收风箬进行微型飞行器（ARMAVs）空中回收的新概念。母船拖曳着连接着电缆的风箬，并控制风箬以匹配MAV的飞行模式。本文使用高斯原理推导出电缆-风箬系统的动态模型。在有风条件下，可控制的风箬在回收MAV中起着关键作用。我们开发了一种利用风箬阻力系数的控制方法。基于多连杆电缆-风箬系统的仿真结果展示了空中回收概念的可行性和风箬的可控性。 

牵引车辆（母船）被命令按照一个圆形轨道行驶，从而使得被牵引物体（风箬）也在一个较小半径和较低速度相对于母船的圆形轨道上运行。要被回收的微型飞行器（MAV）被规划为跟随风箬轨道，并以相对较低的空速接近风箬。在MAV对接风箬后，牵引电缆和物体（风箬和MAV）被缠绕进母船完成空中回收过程。
使用高斯原理推导了电缆-风箬系统的动力学模型。

利用母船和牵引风箬在飞行中回收微型飞行器（MAVs）的方法：基于高斯原理的电缆-风箬系统动力学模型研究

一、研究背景与核心问题

微型飞行器（MAVs）在军事侦察、环境监测等领域具有广泛应用，但其续航能力受限，需通过空中回收技术实现重复使用。传统回收方法（如着陆或网捕）对环境要求高，而利用母船牵引风箬（回收伞）的空中回收技术，通过动态轨道匹配实现低速、低风险对接，成为研究热点。
核心挑战：

1. 非完整约束：风箬通过固定长度的电缆与母船连接，运动受几何约束；
2. 非线性空气动力：风箬受阻力、风扰等复杂力作用，动力学模型高度非线性；
3. 动态耦合：母船、电缆、风箬与MAV的运动相互影响，需协同控制。

二、高斯原理在动力学建模中的应用

高斯原理（最小约束原理）通过虚功原理直接处理约束力，避免传统拉格朗日乘子法的复杂性，适用于非完整约束系统。
建模步骤：

1. 系统简化与假设：
   • 忽略电缆弹性与质量，视为刚性连接；
   • 风箬简化为质点，受空气阻力、重力与电缆张力作用。
2. 坐标系与运动学定义：
   • 以母船拖曳点为原点建立惯性坐标系；
   • 风箬位置向量：r=(x,y,z)，速度：v=(x˙,y˙​,z˙)。
3. 动力学方程推导：

   • 根据高斯原理，系统虚功满足：

其中， 为合外力（阻力、重力、张力）， 为虚位移，需满足约束条件 。

• 引入拉格朗日乘子 λ 处理约束力（电缆张力），得到方程组：

其中， 为阻力系数，为电缆长度， 为风箬质量，为重力加速度。

4. 数值求解方法：

• 采用四阶龙格-库塔法求解非线性微分方程组；
• 关键参数：初始位置 (x0​,y0​,z0​)、速度 (vx0​,vy0​,vz0​)、物理参数 m=10kg,CD​=1.2,l=50m。

三、关键动态特性分析

1. 张力波动与惯性耦合：
   • 母船机动时，电缆张力传递存在时滞效应，转弯时峰值张力可达静载的2倍；
   • 若考虑电缆弹性（杨氏模量 E=9×109N/m2），需引入应变方程 σ=Eϵ，模型扩展为分布参数系统。
2. 空气动力非线性：
   • 风箬面积每增加20%，稳态回收速度下降约15%（仿真验证）；
   • 风扰建模：外部风速 vw​ 引入相对速度项 vrel​=∥v−vw​∥，修正阻力公式。
3. 控制策略集成：
   • 阻力系数调节：通过改变风箬攻角调整 CD​，实现轨道半径控制（如PI控制器使半径收敛至5m内）；
   • 轨迹跟踪：基于李雅普诺夫函数设计控制器，确保MAV稳定跟踪风箬轨道。

四、仿真与验证

1. 典型工况模拟：
   • 母船按圆形轨道飞行，风箬形成较小半径的相对轨道；
   • MAV采用导弹导引策略拦截风箬，仿真显示在 CD​=1.2、初始速度5m/s条件下，MAV可在10秒内稳定对接，位置误差<0.3m。
2. 参数敏感性分析：
   • 电缆直径从0.027m增至0.047m，张力波动降低40%（因刚度提升）；
   • 弹性模量 E 降低50%时，波动幅值增加60%，需阻尼补偿。

五、模型扩展与优化方向

1. 多体耦合建模：

   • 将电缆离散为N段刚体（集中质量法），每段动力学方程为：

其中 为相邻段张力。

2. 高斯过程增强：

• 结合导数和观测数据的高斯过程回归（GPR），提升非线性阻力预测精度（RMSE<5%）。

1. 实时控制集成：

   • 嵌入模型预测控制（MPC），优化目标函数：

其中 为控制输入（如阻力系数、母船速度）。

六、结论与展望

通过高斯原理建立的电缆-风箬系统动力学模型，有效表征了非线性动力学行为，仿真验证了回收概念的可行性。未来工作需进一步融合弹性动力学与鲁棒控制策略，提升复杂环境下的回收可靠性，例如：

• 考虑电缆弹性与风扰的强耦合作用；
• 开发基于深度学习的实时轨迹预测与控制算法。

研究价值：该模型为MAVs空中回收提供了理论支撑，推动了微型飞行器技术的实用化进程。

📚2 运行结果

运行视频：

链接：https://pan.baidu.com/s/1LBLYt8IN9dIXrBSJGyrTHw 
提取码：4cia 
--来自百度网盘超级会员V5的分享

部分代码：

% end mdlInitializeSizes
%
%=======================================================================
% mdlUpdate
% Handle discrete state updates, sample time hits, and major time step
% requirements.
%==============================2=========================================
%
function [xup,fig_joint,fig_cable] = mdlUpdate(t,x,uu,P,fig_joint,fig_cable)

  initialize        = x(1);      % initial graphics
  fig_mothership    = x(2);      % mothership handle
  fig_drogue        = x(3);      % drogue handle
  fig_mav           = x(4);      % mav handle
  
  NN = 0;
  mothership.n     = uu(1+NN);  % north position
  mothership.e     = uu(2+NN);  % east position
  mothership.d     = uu(3+NN);  % down position
  mothership.psi   = uu(4+NN);  % heading angle
  mothership.phi   = uu(5+NN);  % roll angle
  mothership.theta = uu(6+NN);  % flight path angle
  
  NN = 9;  % ignore mothership accelerations
  drogue.n  = uu(1+NN);
  drogue.e  = uu(2+NN);
  drogue.d  = uu(3+NN);
  drogue.vn = uu(4+NN);
  drogue.ve = uu(5+NN);
  drogue.vd = uu(6+NN);
  drogue.phi  = uu(7+NN);
  drogue.theta = uu(8+NN);
  drogue.psi   = uu(9+NN);
  
  NN = NN+9; %
  mav.n     = uu(1+NN);
  mav.e     = uu(2+NN);
  mav.d     = uu(3+NN);
  mav.psi   = uu(4+NN);
  mav.V     = uu(5+NN);
  mav.phi   = uu(6+NN);
  mav.theta = uu(7+NN);
  
  if (P.cable.N>=2)
      % cable states
      for i = 1:(P.cable.N-1)
          NN = 6*(i-1) + 9 + 9 + 9;
          cable.n(i)  = uu(1+NN);
          cable.e(i)  = uu(2+NN);
          cable.d(i)  = uu(3+NN);
          cable.vn(i) = uu(4+NN);
          cable.ve(i) = uu(5+NN);
          cable.vd(i) = uu(6+NN);
      end
      
      X(1,1:2) = [mothership.n, cable.n(1)];
      Y(1,1:2) = [mothership.e, cable.e(1)];
      Z(1,1:2) = [mothership.d, cable.d(1)];

      NN = P.cable.N ;
      X(NN,1:2) = [cable.n(NN-1), drogue.n];
      Y(NN,1:2) = [cable.e(NN-1), drogue.e];
      Z(NN,1:2) = [cable.d(NN-1), drogue.d];
      
      if P.cable.N >= 3
          for i = 2:(P.cable.N-1)
              NN = i ;
              X(NN,1:2) = [cable.n(NN-1), cable.n(NN)];
              Y(NN,1:2) = [cable.e(NN-1), cable.e(NN)];
              Z(NN,1:2) = [cable.d(NN-1), cable.d(NN)];
          end
      end
  end
  
  NN = P.cable.N - 1;
  ell = [...
      cable.n(NN) - drogue.n;...
      cable.e(NN) - drogue.e;...
      cable.d(NN) - drogue.d;...
      ];
  drogue.psi = atan2(ell(2),ell(1));
  
  
    
  if initialize==0, 
      % initialize the plot
      initialize = 1;

      %%%%%%%%%%%%%%%%%%
      % create figure
      figure(1), clf, hold on
      axis([P.x_min, P.x_max, P.y_min, P.y_max, P.z_min, P.z_max]);

      % plot mothership, cable, drogue
      fig_mothership = drawUav(mothership, P.mothership.size, P, [], 'normal');
      fig_drogue     = drawUav(drogue, P.drogue.size, P, [], 'normal');
      fig_mav        = drawUav(mav, P.mav.size, P, [], 'normal');
      
      % Plot cable with multiple links

      for i = 1:P.cable.N
          fig_cable(i) = drawCable(X(i,1:2),Y(i,1:2),Z(i,1:2), P, [], 'normal');
      end
      
      for i = 1:(P.cable.N-1)
          fig_joint(i) = drogueJoint([cable.e(i); cable.n(i); -cable.d(i)], P.joint_size, [], '.red', 'normal');
      end

      grid on
      
      title('Aerial Recovery Demo')
      xlabel('East,m');
      ylabel('North,m');
      zlabel('Altitude,m');
      %view(-154,50)
      view(24,50)

  else  % do this at every time step

      % plot mothership, cable, drogue
      drawUav(mothership, P.mothership.size, P, fig_mothership);
      drawUav(drogue, P.drogue.size, P, fig_drogue);      
      drawUav(mav, P.mav.size, P, fig_mav);
      
%       plot3(mothership.e, mothership.n, -mothership.d);hold on;
%       plot3(drogue.e, drogue.n, -drogue.d);
%       plot3(mav.e, mav.n, -drogue.d);hold on;      
      
      for i = 1:P.cable.N
          drawCable(X(i,1:2),Y(i,1:2),Z(i,1:2), P, fig_cable(i));
      end
      
      for i = 1:(P.cable.N-1)
          drogueJoint([cable.e(i); cable.n(i); -cable.d(i)], P.joint_size, fig_joint(i));
      end
      drawnow
%      drawBreadcrumb(mothership.n,mothership.e,mothership.d,'r');
%      drawBreadcrumb(drogue.n,drogue.e,drogue.d,'g');
      
      % plot the trails of the system
      plot3(mothership.e, mothership.n, -mothership.d);hold on;
      plot3(drogue.e, drogue.n, -drogue.d);hold on;
      plot3(mav.e, mav.n, -mav.d);hold on; 

  end
  xup(1) = initialize;
  xup(2) = fig_mothership;
  xup(3) = fig_drogue;
  xup(4) = fig_mav;
%   xup = [initialize; fig_mothership; fig_drogue; fig_mav];

%end mdlUpdate

%=======================================================================
% mdlOutput
%==============================2=========================================
%
function sys = mdlOutput(t,x,u,P)

  sys = [];
  
%----------------------------------------------------------------------
%----------------------------------------------------------------------
% User defined functions
%----------------------------------------------------------------------
%----------------------------------------------------------------------

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function handle =drawCable(X, Y, Z, P, handle, mode);
  % cablePlot:  plot cable between mothership and drogue
  
  if isempty(handle)
      handle = plot3(Y,X,-Z,...
          'color','b',...
          'EraseMode', mode);
  else
      set(handle,'XData',Y,'YData',X,'ZData',-Z);
%       drawnow
  end

  
function handle = drogueJoint(z, size, handle, color, mode)

  if isempty(handle),
    handle = plot3(z(1), z(2), z(3), color, 'EraseMode', mode);
  else
    set(handle,'XData',z(1),'YData',z(2),'ZData',z(3));
%     drawnow
  end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function handle=drawUav(uav, size, P, handle, mode);
  % uavPlot:  plot UAV attitude at Euler angles phi, theta, psi

  [Vert_uav,Face_uav,colors_uav] = uavVertFace(P);
  % rescale vertices
  Vert_uav = size*Vert_uav;
  % rotate vertices by phi, theta, psi
  Vert_uav = rotateVert(Vert_uav, pi-uav.phi, -uav.theta, uav.psi);
  % transform vertices from NED to XYZ
  R = [...
      0, 1, 0;...
      1, 0, 0;...
      0, 0, 1;...
      ];
  Vert_uav = Vert_uav*R;
  % translate vertices in XYZ
  Vert_uav = translateVert(Vert_uav, [uav.e; uav.n; -uav.d]);

  % collect all vertices and faces
  V = [...
      Vert_uav;...
      ];
  F = [...
      Face_uav;...
      ]; 
  patchcolors = [...
      colors_uav;...
      ];

  if isempty(handle),
  handle = patch('Vertices', V, 'Faces', F,...
                 'FaceVertexCData',patchcolors,...
                 'FaceColor','flat',...
                 'EraseMode', mode);
  else
    set(handle,'Vertices',V,'Faces',F);
%     drawnow
  end

 
%%%%%%%%%%%%%%%%%%%%%%%
function [V,F,patchcolors] = uavVertFace(P);
  % uavData:  define patch faces of uav

 
  % vertices of the fuselage
  Vert_fuse = [...
        P.drawuav.l_fuse/2, P.drawuav.w_fuse/2, P.drawuav.h_fuse/2;...
        P.drawuav.l_fuse/2, P.drawuav.w_fuse/2, -P.drawuav.h_fuse/2;...
        P.drawuav.l_fuse/2, -P.drawuav.w_fuse/2, -P.drawuav.h_fuse/2;...
        P.drawuav.l_fuse/2, -P.drawuav.w_fuse/2, P.drawuav.h_fuse/2;...
        -P.drawuav.l_fuse/2, P.drawuav.w_fuse/2, -P.drawuav.h_fuse/2;...
        -P.drawuav.l_fuse/2, P.drawuav.w_fuse/2, P.drawuav.h_fuse/2;...
        -P.drawuav.l_fuse/2, -P.drawuav.w_fuse/2, -P.drawuav.h_fuse/2;...
        -P.drawuav.l_fuse/2, -P.drawuav.w_fuse/2, P.drawuav.h_fuse/2];    
  % define faces of fuselage
  Face_fuse = [...
        1, 2, 3, 4;... % front
        5, 6, 8, 7;... % back
        2, 3, 7, 5;... % top
        1, 2, 5, 6;... % right 
        3, 4, 8, 7;... % left
        1, 4, 8, 6];   % bottom

  Vert_wing_r = [...
        P.drawuav.l_fuse/2-P.drawuav.w_recess, P.drawuav.w_fuse/2, P.drawuav.w_height/2;...
        P.drawuav.l_fuse/2-P.drawuav.w_recess, P.drawuav.w_fuse/2, -P.drawuav.w_height/2;...
        -P.drawuav.l_fuse/2, P.drawuav.w_fuse/2, 0;...
        -P.drawuav.w_back, P.drawuav.w_fuse/2+P.drawuav.w_width, P.drawuav.w_height/2;...
        -P.drawuav.w_back, P.drawuav.w_fuse/2+P.drawuav.w_width, -P.drawuav.w_height/2;...
        -P.drawuav.w_back-P.drawuav.w_tipwidth, P.drawuav.w_fuse/2+P.drawuav.w_width, 0;...
    ];

  Face_wing_r = 8 + [...
        1, 2, 5, 4;... % front
        4, 5, 6, 4;... % side
        2, 3, 6, 5;... % top
        1, 3, 6, 4;... % bottom
    ];

  Vert_wing_l = [...
        P.drawuav.l_fuse/2-P.drawuav.w_recess, -P.drawuav.w_fuse/2, P.drawuav.w_height/2;...
        P.drawuav.l_fuse/2-P.drawuav.w_recess, -P.drawuav.w_fuse/2, -P.drawuav.w_height/2;...
        -P.drawuav.l_fuse/2, -P.drawuav.w_fuse/2, 0;...
        -P.drawuav.w_back, -P.drawuav.w_fuse/2-P.drawuav.w_width, P.drawuav.w_height/2;...
        -P.drawuav.w_back, -P.drawuav.w_fuse/2-P.drawuav.w_width, -P.drawuav.w_height/2;...
        -P.drawuav.w_back-P.drawuav.w_tipwidth, -P.drawuav.w_fuse/2-P.drawuav.w_width, 0;...
    ];

  Face_wing_l = 14 + [...
        1, 2, 5, 4;... % front
        4, 5, 6, 4;... % side
        2, 3, 6, 5;... % top
        1, 3, 6, 4;... % bottom
    ];

  Vert_rud_r = [...
        -P.drawuav.w_back, P.drawuav.w_fuse/2+P.drawuav.w_width, P.drawuav.w_height/2;...
        -P.drawuav.w_back-P.drawuav.w_tipwidth, P.drawuav.w_fuse/2+P.drawuav.w_width, 0;...
        -P.drawuav.w_back-P.drawuav.w_tipwidth, P.drawuav.w_fuse/2+P.drawuav.w_width, -P.drawuav.r_height;...
        -P.drawuav.w_back-P.drawuav.w_tipwidth+P.drawuav.r_top, P.drawuav.w_fuse/2+P.drawuav.w_width, -P.drawuav.r_height;...
    ];

  Face_rud_r = 20 + [...
        1, 2, 3, 4;...
    ];

  Vert_rud_l = [...
        -P.drawuav.w_back,              -P.drawuav.w_fuse/2-P.drawuav.w_width, P.drawuav.w_height/2;...
        -P.drawuav.w_back-P.drawuav.w_tipwidth, -P.drawuav.w_fuse/2-P.drawuav.w_width, 0;...
        -P.drawuav.w_back-P.drawuav.w_tipwidth, -P.drawuav.w_fuse/2-P.drawuav.w_width, -P.drawuav.r_height;...
        -P.drawuav.w_back-P.drawuav.w_tipwidth+P.drawuav.r_top, -P.drawuav.w_fuse/2-P.drawuav.w_width, -P.drawuav.r_height;...
    ];

  Face_rud_l = 24 + [...
        1, 2, 3, 4;...
    ];

V = [Vert_fuse; Vert_wing_r; Vert_wing_l; Vert_rud_r; Vert_rud_l];
F = [Face_fuse; Face_wing_r; Face_wing_l; Face_rud_r; Face_rud_l];

patchcolors = [...
    P.myblue;... % fuselage front
    P.myblue;... % back
    P.myyellow;... % top
    P.myblue;... % right
    P.myblue;... % left
    P.myblue;... % bottom
    P.myblue;... % right wing front
    P.mygreen;... % side
    P.myblue;...  % top
    P.myred;...   % bottom
    P.myblue;... % left wingfront
    P.mygreen;... % side
    P.myblue;...  % top
    P.myred;...   % bottom
    P.mygreen;...% right rudder
    P.mygreen;...% left rudder
    ];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cityPlot(P)

    for i=1:P.num_blocks,
      for j=1:P.num_blocks,
        C = [...
            i*P.street_width + (i-1)*P.building_width + 0.5*P.building_width;...
            j*P.street_width + (j-1)*P.building_width + 0.5*P.building_width;...
            ];
        X = [...
            C(1) - 0.5*P.building_width;...
            C(1) - 0.5*P.building_width;...
            C(1) + 0.5*P.building_width;...
            C(1) + 0.5*P.building_width;...
            ];
        Y = [...
            C(2) - 0.5*P.building_width;...
            C(2) + 0.5*P.building_width;...
            C(2) + 0.5*P.building_width;...
            C(2) - 0.5*P.building_width;...
            ];
        fill(X,Y,'g')
      end
    end
    
%%%%%%%%%%%%%%%%%%%%%%%
function Vert=rotateVert(Vert,phi,theta,psi);
% Rotation matrix from body coordinates to vehicle coordinates
 

🎉3 参考文献

文章中一些内容引自网络，会注明出处或引用为参考文献，难免有未尽之处，如有不妥，请随时联系删除。

 

🌈4 Matlab代码、数据、文章

资料获取，更多粉丝福利，MATLAB|Simulink|Python资源获取