基于PID调节的简易无人机飞行高度控制系统分析

Posted on In

  自动控制原理实验报告

基于PID调节的简易无人机飞行高度控制系统分析

实验背景

  无人机技术在航拍测绘、物流配送、农业植保、应急救援等领域的普及，使其飞控系统的可靠性成为技术落地的核心前提 —— 而高度控制是飞控系统的 “基础生命线”：若无人机无法稳定维持目标高度，不仅会影响航拍画质、物流投送精度等功能实现，甚至会因高度失控引发坠机风险。

  实际工业级飞控系统需融合姿态解算、多传感器融合、抗干扰算法等复杂模块，直接上手开发门槛较高。因此，本实验采用垂直方向动力学简化模型（忽略姿态耦合、风场强扰动等复杂因素），聚焦 “高度控制” 这一核心子问题，是入门级飞控算法验证的典型载体。

  传统无PID控制器下，无人机高度系统易出现 “响应迟缓、超调剧烈、稳态误差大” 的问题，无法满足基本悬停需求。而 PID 控制器作为工业控制领域最经典、易实现的闭环控制算法，兼具 “结构简单、参数物理意义明确、适配性强” 的优势，是验证高度控制算法有效性的理想选择。

  本实验通过 “简化模型 + PID 控制” 的组合，既贴合入门级飞控算法的学习逻辑，也能直观验证闭环控制对系统性能的优化价值，为后续复杂飞控系统的设计与调试奠定理论与实践基础。

数学建模

  忽略无人机的飞行姿态，仅考虑垂直方向动力学：

  在垂直方向上无人机受到升力\(F\)、自身重力\(G\)和空气阻力\(f\)。受力情况如图1所示。


图1：无人机受力情况

  在考虑增量情况下进行受力分析得到：

\[ \Delta F - kv = ma \]

  其中\(k\)为空气阻力系数、\(\Delta F\)为推力增量。将速度和加速度代入上式得到：

\[ \Delta F(t) - k \frac{dx(t)}{dt} = m \frac{d^2x(t)}{dt^2} \]

  无人机质量为\(m\)，空气阻力系数为\(k\)。在零初始条件下，对上述方程各项进行拉普拉斯变换得到复数域数学模型：

\[ ms^2X_0(s)+ksX_0(s) = \Delta F'(s) \]

  由定义得到系统的开环传递函数：

\[ G(s) = \frac{X_0(s)}{\Delta F'(s)} = \frac{1}{ms^2+ks} \]

  闭环传递函数:

\[ \Phi(s) = \frac{G(s)}{1+G(s)} = \frac{\frac{1}{ms^2+ks}}{1+\frac{1}{ms^2+ks}} = \frac{1}{ms^2+ks+1} \]

  将闭环传递函数化为二阶系统的标准形式:

\[ \Phi(s) = \frac{\frac{1}{m}}{s^2 + \frac{k}{m}s+\frac{1}{m}} \]

  其中无阻尼振荡频率\(\omega_n=\sqrt{\frac{1}{m}}\)，阻尼比\(\zeta=\frac{k}{2\sqrt{m}}\)。由于空气阻力系数\(k\)较小，因此可以将阻尼比\(\zeta\)的值视为\(0\)到\(1\)之间，所以该系统为二阶欠阻尼系统。

控制目标

  1. 实现高度阶跃的跟踪，从\(0\rightarrow1m\)悬停

  2. 动态性能指标：

    上升时间短：\(t_r<1s\)

    调节时间短：\(t_s<3s\)

    超调量小：\(\sigma\%<5\%\)

    稳态误差小

  3. 系统稳定，无发散或者持续振荡

无PID控制下的时域分析

  假设无人机质量\(m=5kg\)，空气阻力系数\(k=0.5\)。在这种情况下使用MATLAB的Simulink绘制简单的闭环系统阶跃响应图像。


图2：无PID调节的Simulink控制框图


图3：无PID调节的单位阶跃响应

动态性能指标计算

  根据图2模拟出来的数据可以计算其动态性能指标。对于无人机的飞行高度控制，我们更加关注无人机能否快速地准确地到达我们指定的高度。因此，我们需要计算该系统的上升时间\(t_r\)、超调量\(\sigma\%\)以及调节时间\(t_s\)。

\[ 上升时间：t_r=\frac{\pi-\beta}{\omega_d}=\frac{\pi-\arccos\zeta}{\omega_n\sqrt{1-\zeta^2}} \]

\[ 超调量：\sigma\%=e^{-\pi\zeta/\sqrt{1-\zeta^2}}\times100\% \]

\[ 调节时间：t_s=\frac{3.5}{\zeta\omega_n} ~~~~~~~~(误差带\Delta=0.05) \]

  将无阻尼振荡频率\(\omega_n=\sqrt{\frac{1}{m}}=\sqrt{\frac{1}{5}}\approx0.4472\)和阻尼比\(\zeta=\frac{k}{2\sqrt{m}}=\frac{0.5}{2\sqrt{5}}\approx0.1118\)代入上式得到：

\[ t_r\approx7.069s,~\sigma\%\approx70.23\%,~t_s\approx70.00s \]

  从图3以及动态性能指标的计算结果中可以看出在无PID调节的情况下，虽然最后无人机的高度会趋向于1m，但是响应速度较慢，调节时间较长，并且超调量很高，振荡幅度很大。

PID控制器设计

  这并不是我们希望得到的效果。我们希望无人机在接收到一个0到1m的阶跃信号之后，能够迅速飞行至1m并保持误差可接受的稳定状态。因此需要设计PID控制器改善系统的动态性能。

PID原理

  PID控制的核心逻辑是根据“期望输出和实际输出的误差”，通过比例、积分、微分三个环节的加权组合生成控制器，使得系统能够快速、稳定地跟踪目标值。

  在本实验中我们需要让无人机快速、无超调地跟踪期望高度。而实际高度与期望高度存在误差，需要通过调整推力增量来消除误差。这里PID控制器地作用就是将高度误差转化为推力增量指令，通过P、I、D三个环节的协同作用来快速消除误差避免系统振荡。

比例环节P

  在比例环节中，其输出与系统的当前误差成正比。我们记比例环节的输出为\(u_P\)，比例系数为\(K_p\)，实时误差为\(e(t)\)，则有如下式子成立：

\[ u_P=K_p\cdot e(t) \]

  其中，系统的实时误差由期望高度减去实际高度得到。以上式子表示的物理含义为：当前的误差有多大，就给出多大的控制量。也就是说误差越大，无人机的推力就越强烈。

  若无人机的实际高度低于1m，则有\(e(t)>0\)，P环节就会输出正的推力增量，电机加大推力，无人机上升；若无人机的实际高度高于1m，则有\(e(t)<0\)，P环节就会输出负的推力增量，无人机就会下降。

积分环节I

  在积分环节中，其输出与累计的误差值成正比。它反映了误差的历史总和。我们记积分环节的输出为\(u_I\)，积分系数为\(K_i\)，则有以下式子成立：

\[ u_I = K_i \cdot \int_0^t e(\tau)d\tau \]

  以上式子表示的物理含义为：只要误差存在，就持续累积控制量，直到误差为0。积分环节主要针对于P环节无法消除的稳态误差。

微分环节D

  在微分环节中，其输出与误差的变化率成正比。它反映了误差的变化趋势。我们记微分环节输出为\(u_D\)，微分系数为\(K_d\)，则有一下式子成立：

\[ u_D=K_d\cdot\frac{de(t)}{dt} \]

  以上式子表示的物理含义为：预判误差的变化趋势，提前给出反向控制量，它能够抑制系统超调，让响应更加平稳。当无人机在快速上升过程中，误差\(e(t)\)快速减小，\(\frac{de(t)}{dt} < 0\)，则微分环节输出负的推力增量，提前减小电机推力，避免无人机过冲。

PID控制器的整体模型

  PID的总输出是比例、积分、微分三个环节的叠加，即控制量由以下式子得到：

\[ u(t)=u_P+u_I+u_D=K_p\cdot e(t)+K_i\cdot \int_0^te(\tau)d\tau+K_d\cdot\frac{de(t)}{dt} \]

  在零初始条件下对上式做拉普拉斯变换得到：

\[ U(s) = K_pE(s) + \frac{K_i}{s}E(s) + K_d\cdot sE(s) \]

\[ G_{PID}(s)=\frac{U(s)}{E(s)}=K_p+\frac{K_i}{s}+K_d\cdot s=\frac{K_ds^2+K_ps+K_i}{s} \]

  其中\(G_{PID}(s)\)为PID环节的传递函数。结合之前的无人机高度控制系统的开环传递函数\(G(s)=\frac{1}{ms^2+ks}\)，得到系统的开环传递函数为：

\[ G_{open}(s)=G_{PID}(s)G(s)=\frac{K_ds^2+K_ps+K_i}{ms^3+ks^2}=\frac{K_ds^2+K_ps+K_i}{5s^3+0.5s^2} \]

  闭环传递函数为：

\[ \Phi_{PID}(s)=\frac{G_{PID}(s)\cdot G(s)}{1+G_{PID}(s)\cdot G(s)}=\frac{K_ds^2+K_ps+K_i}{ms^3+(k+K_d)s^2+K_ps+K_i}=\frac{K_ds^2+K_ps+K_i}{5s^3+(0.5+K_d)s^2+K_ps+K_i} \]

  在MATLAB的Simulink中绘制PID控制器框图如下：


图4：含有PID调节的Simulink控制框图

根轨迹分析：PID参数调整

  根轨迹分析能够将闭环极点随PID参数变化的规律进行可视化，通过主动调整\(K_p\)、\(K_i\)、\(K_d\)，将闭环极点配置到s平面的理想区域，从而达到精准控制超调、调节时间等时域指标。

  参数的调整根据\(P\rightarrow PD \rightarrow PID\)的递进顺序，每一步使用根轨迹确定参数范围。

仅比例控制：\(K_p\)参数调节

  在仅比例控制的情况下，令\(K_i=0\)、\(K_d=0\)，系统开环传递函数简化为：

\[ G_{open,P}(s)=\frac{K_p}{s(5s+0.5)}=\frac{\frac{K_p}{5}}{s(s+0.1)} \]

  其中开环系统根轨迹增益\(K^*=\frac{K_p}{5}\)。系统闭环传递函数简化为：

\[ \Phi_P(s)=\frac{K_p}{5s^2+0.5s+K_p} \]

  于是，特征方程式可写为：

\[ 5s^2+0.5s+K_p = 0 \]

  得到特征方程式的根为：

\[ s_1=-0.05+0.1\sqrt{0.25-20K_p} \]

\[ s_2=-0.05-0.1\sqrt{0.25-20K_p} \]

  在MATLAB中可以绘制该系统当\(K_p\)从零到无穷时候的根轨迹。编写如下代码：

num = 1;
den = conv([1, 0], [5, 0.5]);
sys = tf(num, den);
rlocus(sys);

  得到根轨迹图像如图:


图5：比例控制下的根轨迹图

  根据上图分析可以发现，当增益\(K_p<0.0125\)时候，闭环特征方程两根分布在实轴上，系统为过阻尼系统，无超调。虽然超调量符合系统设计要求，但是在这种情况下系统响应较慢。

  而上升时间\(t_r\)与阻尼角\(\beta\)和有阻尼振荡频率\(\omega_d\)有关，根据公式\(t_r=\frac{\pi-\beta}{\omega_d}\)可知：为使系统响应尽可能快，则上升时间要尽可能小，所以阻尼角和无阻尼振荡频率要尽可能大。根据图6所示的欠阻尼二阶系统的特征参量可以很容易发现，我们需要让闭环特征根虚部的绝对值尽可能大。


图6：欠阻尼二阶系统的特征参量

  所以根据以上分析可以得到，在仅使用比例控制的情况下，要使得系统能够快速响应，则参数\(K_p\)值要适当大一些。目前我们选择\(K_p=20\)。


图7：不同Kp对响应的影响

比例-微分控制：\(K_d\)参数调节

  在\(K_p=20\)的情况下，虽然系统的上升时间\(t_r\)变短了，但是超调量较大并且调节时间\(t_s\)较长。因此在比例控制的基础上，需要再使用微分控制来抑制系统超调，让响应更加平稳。

  当\(K_p=20,~K_i=0\)时，系统的开环传递函数为：

\[ G_{open,PD}(s) = \frac{K_ds^2+20s}{5s^3+0.5s^2}=\frac{K_ds+20}{5s^2+0.5s} \]

  现在需要对其绘制根轨迹以便确定\(K_d\)的范围。但是由于系统开环传函中\(K_d\)是s项的系数，而非整个开环传函的增益。因此不能在MATLAB中使用常规的 rlocus()函数来绘制根轨迹。

  为了绘制根轨迹，需要参数化遍历\(K_d\)，逐个计算闭环特征方程的极点并在复平面上绘制。需要编写以下代码：

acc = 100;  % 精度
Kd_range = linspace(0, 30, acc);  % Kd范围从0~30
poles_all = zeros(length(Kd_range), 2);  % 二阶系统，每个K_d对应2个闭环极点

% 遍历K_d，计算每个K_d对应的闭环极点
for i = 1:length(Kd_range)
    Kd = Kd_range(i);
    % 构造开环传函和单位负反馈闭环系统
    num = [Kd 20];    % 分子：Kd*s + 20
    den = [5 0.5 0];  % 分母：5s² + 0.5s
    G_open = tf(num, den);
    sys_cl = feedback(G_open, 1);
    poles_all(i, :) = pole(sys_cl);
end

% 绘制根轨迹
plot(real(poles_all), imag(poles_all), 'b-', 'LineWidth',1.2);  % 极点轨迹
hold on;
plot(real(poles_all), -imag(poles_all), 'b-', 'LineWidth',1.2); % 共轭极点

% 辅助标注
grid on;
xlabel('实轴'); ylabel('虚轴');
title('G_{open,PD}(s)的根轨迹');
plot(real(poles_all(1,:)), imag(poles_all(1,:)), 'ro', 'MarkerSize',6); % K_d=0
text(real(poles_all(1,1))+0.1, imag(poles_all(1,1)), 'K_d=0');
plot(real(poles_all(end,:)), imag(poles_all(end,:)), 'go', 'MarkerSize',6); % K_d=30
text(real(poles_all(end,1))+0.1, imag(poles_all(end,1)), 'K_d=30');

  运行得到如下图像：


图8：精度为100下的根轨迹

  观察上图会发现在分离点附近出现一个菱形，实际上这是\(K_d\)精度不足导致的。精度为\(100\)时候步长较大，在分离点附近跨过了分离点的值。为了根轨迹看上去更加美观，可以尝试提高精度。当精度 acc设置为\(10000\)时候，几乎就看不出菱形了。


图9：精度为10000下的根轨迹


图10：根轨迹分离点附近的Kd值

  观察图像可以发现，当\(K_d<19.5\)时，闭环特征方程有两个复数根，系统为欠阻尼，存在振荡。当\(K_d>19.5\)时，闭环特征方程有两个实数根，系统为过阻尼，无振荡。

  为了找到合适的\(K_d\)参数，需要根据系统闭环传递函数的式子以及动态时域性能指标进行分析。在PD控制的情况下，系统的闭环传递函数为：

\[ \Phi_{PD}(s)=\frac{K_ds+20}{5s^2+(0.5+K_d)s+20} \]

可以发现，该系统为带零点的二阶系统，其手动计算性能指标较为复杂。因此可以使用MATLAB中的 step()和 stepinfo()函数直接进行仿真计算。编写如下代码：

%构造闭环传递函数
num = [Kd 20];        % 分子：K_d s + 20
den = [5 0.5+Kd 20];  % 分母：5s² + (0.5+K_d)s + 20
sys = tf(num, den);

step(sys);  % 单位阶跃响应
hold on;
grid on;
xlabel('时间/s');
ylabel('高度/m');
title('PD控制系统单位阶跃响应');

info = stepinfo(sys);  % 提取阶跃响应性能指标

% 输出结果
fprintf('===== 性能指标（K_d=%.2f） =====\n', Kd);
fprintf('上升时间：%.2f s（10%%→90%%稳态值）\n', info.RiseTime);
fprintf('超调量：%.2f %%\n', info.Overshoot);
fprintf('调节时间：%.2f s（±5%%误差带）\n', info.SettlingTime);

  运行代码，分别绘制当\(K_d=5,K_d=10,K_d=19.5,K_d=25\)时系统的单位阶跃响应。


图11：不同Kd值的单位阶跃响应


图12：不同Kd值单位阶跃响应的时域性能

  可以发现，PD控制中微分环节对比例环节的振荡效果有明显的抑制作用。并且随着\(K_d\)的增大，抑制效果越强，同时系统超调量、上升时间和调节时间均减少，系统的响应速度更快。为了满足\(t_r<1s,~t_s<3s,~\sigma\%<5\%\)的控制目标，可以选择\(K_d=35\)。


图13：Kd=35时系统时域性能指标

比例-积分-微分控制：\(K_i\)参数调节

  从上面的仿真结果来看，PD控制已经基本符合条件。而在PID中加入积分环节的作用是为了消除稳态误差。我们可以尝试分析一下在\(K_p,~K_d\)两个参数固定时候，\(K_i\)对系统稳定性的影响。

  当\(K_p=20,~K_d=35\)时，系统的闭环传函为：

\[ \Phi_{PID}(s)=\frac{35s^2+20s+K_i}{5s^3+35.5s^2+20s+K_i} \]

  这是一个高阶系统，手动计算分析难度较大，因此可以借助MATLAB来绘制根轨迹以及随着\(K_i\)的变化其阶跃响应的变化。编写以下代码来绘制根轨迹：

Ki_range = linspace(-1, 300, 100000);  % K_i从-1到300
poles_all = zeros(length(Ki_range), 3);

for i = 1:length(Ki_range)
    Ki = Ki_range(i);
    coeffs = [5, 35.5, 20, Ki];  % 闭环特征方程：5s^3 + 35.5s^2 + 20s + Ki = 0
    poles = roots(coeffs);     
    poles_all(i, :) = poles;  
end

hold on; grid on; axis equal;

% 绘制3条根轨迹分支
plot(real(poles_all(:,1)), imag(poles_all(:,1)), 'b-', 'LineWidth',1.2);
plot(real(poles_all(:,2)), imag(poles_all(:,2)), 'r-', 'LineWidth',1.2);
plot(real(poles_all(:,3)), imag(poles_all(:,3)), 'g-', 'LineWidth',1.2);
plot([0,0], [-50,50], 'k--', 'LineWidth',1);

% 标注K_i=0时的初始极点
poles_Ki0 = roots([5,35.5,20,0]);
plot(real(poles_Ki0), imag(poles_Ki0), 'ko', 'MarkerSize',8);
text(real(poles_Ki0), imag(poles_Ki0), 'K_i=0', 'FontSize',9);

% 标题与轴标签
xlabel('实轴'); ylabel('虚轴');
title('PID系统根轨迹');
legend('极点分支1','极点分支2','极点分支3','Location','best');


图14：PID系统根轨迹


图15：临界点处的Ki值

  运行得到如上图所示的图像。可以发现，随着\(K_i\)的增大，系统并不会变得更加稳定，反而会变得更加不稳定。当\(K_i>141.856\)时，系统的闭环特征方程根会出现在右半平面上，此时系统是发散的。


图16：不同Ki的单位阶跃响应

  结合图16以及以上分析，\(K_i\)的值不能取太大。为满足控制目标，可以取\(K_i=1\)，甚至取\(K_i=0\)去掉积分环节都可以。


图17：Ki为1时的时域性能指标

频域分析：Nyquist图与Bode图绘制

  现在我们已经确定了PID三个参数的值：

\[ K_p=20,~~K_i=1,~~K_d=35 \]

  此时系统开环传函为：

\[ G(s)=\frac{35s^2+20s+1}{5s^3+0.5s^2} \]

  闭环传函为：

\[ \Phi_{PID}(s)=\frac{35s^2+20s+1}{5s^3+35.5s^2+20s+1} \]

  可以根据开环传函分析其频率特性。将\(s=j\omega\)代入得到：

\[ G(j\omega)=\frac{35\omega^2-1-20\omega j}{0.5\omega^2+5\omega^3j}=\frac{\sqrt{(35\omega^2-1)^2+(20\omega)^2}e^{j\arctan{\frac{20\omega}{1-35\omega^2}}}}{\sqrt{(0.5\omega^2)^2+(5\omega^3)^2}e^{j\arctan{10\omega}}} \]

  所以系统的幅频特性和相频特性分别为：

\[ A(\omega)=\frac{1}{5\omega^2}\sqrt{\frac{(35\omega^2-1)^2+400\omega^2}{0.01+\omega^2}} \]

\[ \phi(\omega)=-\pi + \begin{cases} \arctan{\frac{20\omega}{1-35\omega^2}}-\arctan10\omega & 0<\omega<\frac{1}{\sqrt{35}} \\ \pi-\arctan{\frac{20\omega}{35\omega^2-1}}-\arctan10\omega & \omega>\frac{1}{\sqrt{35}} \end{cases} \]

Nyquist图绘制

  基于系统的频率特性，可以很容易在MATLAB中画出其奈氏图。运行如下代码：

num = [35 20 1];
den = [5 0.5 0 0];
G = tf(num, den);

nyquistplot(G); 
hold on;
plot(-1, 0, 'rx', 'MarkerSize', 10, 'LineWidth', 2);
text(-1, 0.5, '(-1,j0)临界点', 'FontSize', 10);
title('开环传函奈奎斯特图');
xlabel('实轴'); ylabel('虚轴');

  得到系统的奈氏图。


图18：奈奎斯特图

  可以根据上图结合奈奎斯特判据来验证系统的稳定性。正负穿越表示的奈奎斯特判据的表达式为：

\[ Z=P+2(N_--N_+) \]

根据奈奎斯特图像可以看出\(N_-=0,~N_+=0\)，并且根据之前的分析我们知道\(P=0\)。因此\(Z=0\)，所以系统是稳定的。

Bode图绘制

  在MATLAB中使用类似的方法，可以绘制系统的Bode图。运行以下代码：

num = [35 20 1];
den = [5 0.5 0 0];
G = tf(num, den);

bode(G);
grid on; 
hold on;
title('开环传递函数Bode图');
xlabel('频率 (rad/s)');

  得到系统的Bode图


图19：伯德图

  结合先前计算得到幅频特性表达式\(A(\omega)\)和相频特性表达式\(\phi(\omega)\)可以看出计算没有问题。

优化程度分析与总结

  对于这个实验假设的简易无人机飞行高度控制系统，我们使用了PID控制器让系统对于单位阶跃响应能够快速达到控制目标。其中参数分别为：

\[K_p=20,~~K_i=1,~~K_d=35\]

  现在可以结合无PID情况对比一下加上PID控制器以后的系统时域性能指标。

性能指标 无PID 有PID 优化程度
上升时间\(t_r\) 7.069s 0.27s 快了26倍
超调量\(\sigma\%\) 70.23% 4.96% 小了14倍
调节时间\(t_s\) 70s 2.67s 快了26倍


图20：有无PID的对比

  结合表格以及图20中的数据，无 PID 的开环系统呈现强振荡、慢收敛的不稳定特性，完全无法满足无人机高度悬停的工程需求；而 PID 通过 “比例（快速响应误差）+ 积分（消除稳态偏差）+ 微分（抑制超调）” 的协同作用，将系统从 “不可用” 优化为 “稳定可靠”，验证了 PID 控制器在低空无人机飞控中的核心适配性。

  本次实验以 “无人机高度控制” 为载体，明确了 PID 控制器的优化价值：在简易被控对象下，通过整定参数。可同时实现 “快速响应、低超调、短调节时间” 的时域性能，为实际无人机飞控系统的控制器设计提供了基础验证与参数整定思路。