推广

[控制系列] 二、自动控制的本质

iseeyu2年前 (2024-08-24)推广103

我们最常用到的控制理论是pid控制，但是仅仅掌握经典控制理论可能是不够的。因为，现在机器人系统的复杂程度和算力在不断提高，更多先进的控制理论的使用也越来越广泛，所以要不断去深入学习像现代控制理论、非线性控制理论、强化学习等。像一些无人驾驶、imu融合、四足机器人，会涉及到现代控制理论，优化控制，LQR，MPC。；类似人形机器人的这类复杂系统，还可能会涉及到非线性控制。

但是经典控制理论依然是基础，也有广泛应用的基础，所以我们首先要能够充分掌握。

接下来我们一起来学习经典控制理论，同时我会尽可能多的穿插仿真，最后我也会大家一起搭建简单的机器人，一起来感受一下控制的魅力。

在开始之前，我们先来用一个简单的例子看看什么是控制系统建模，控制模型在系统中究竟都在干什么，怎么达到我们期望的效果呢？

2.1 无阻尼钟摆运动

我们以典型的钟摆系统为例，如下为无阻尼钟摆系统

无阻尼单摆系统

根据牛顿第二定律，给出系统的力平衡方程，切向的力与切向加速度，由如下关系

mL^2\ddot \theta + mgL sin\theta = 0\\

等式两边同时除以 mL^2 进行简化得到下面式子：

\ddot \theta+\frac{g}{L}sin\theta=0\\

式中， \theta 为钟摆偏离铅垂线的角度，g是重力加速度，L是摆长。选择变量 x_1=\theta,x_2=\dot \theta ，这样一个二阶微分方程，化简成了一个一阶微分方程组： {\dot{x}}_1=x_2=\dot \theta \\ {\dot{x}}_2=\ddot \theta=-\frac{g}{L}sinx_1\\ 这种非线性钟摆方程，难以利用解析法求解，通过数值求解可以得到一条解的轨迹（这和初始点选择有关），给定一个初始起点，通过matlab数值计算可以得到一条状态轨迹，如下图蓝色的线，其中横坐标 x_1 表示摆的摆角，纵坐标 x_2 为角速度（可想而知，当初始点分布在空间的任何一个位置，那么其解也可能会填满整个空间）

解的一条轨迹

经常将坐标轴（ x_1,x_2 ）称为位置、速度，统称为状态，将由状态组成的空间称为状态空间。从这个图，我们可以看从开始的位置 x_1 趋向于0，速度 x_2 开始增大；在起点与纵轴交点处位置为处0，速度最大；接着位置趋向于最右端，速度开始减小，位置逐渐增大。

注：matlab代码轨迹求解
% 时间范围初值
[t,x]=ode23(@pendulum,[0 5],[-2 -0.25]);
figure%相图解的轨迹
plot(x(:,1),x(:,2))
hold on
plot([min(x(:,1)),max(x(:,1))]*1.2,[0,0],'k.-')
plot([0,0],[min(x(:,2)),max(x(:,2))]*1.2,'k.-')
axis equal
%微分方程无阻尼
function dxdt=pendulum(t,x)
dxdt=[x(2);-sin(x(1))];
end

2.2 有阻尼钟摆运动

从上图，无阻尼钟摆的轨迹可以看出，左右两边的轨迹是对称的，摆会无限往复运动动；当给系统增加阻尼 b ，会怎么样呢！方程如下 {\dot{x}}_1=x_2\\ {\dot{x}}_2=-\frac{g}{L}sinx_1-bx_2\\ 给定一个初始点，数值计算出一条解的轨迹，如下图

有阻尼的解轨迹

注：matlab代码数值计算轨迹
% 时间范围初值
[t,x]=ode23(@pendulum,[0 15],[-2 -0.25]);
figure%相图解的轨迹
plot(x(:,1),x(:,2))
hold on
plot([min(x(:,1)),max(x(:,1))]*1.2,[0,0],'k.-')
plot([0,0],[min(x(:,2)),max(x(:,2))]*1.2,'k.-')
axis equal
%微分方程有阻尼
function dxdt=pendulum(t,x)
dxdt=[x(2);-sin(x(1))-0.5*x(2)];
end

可以看到，轨迹不再左右对称，而是慢慢趋向于零点，摆来回摆动，高度会越来越低，最后趋向静止于最低点。从能量的角度来看，无阻尼系统中没有消耗能量的部分，系统动能、重力势能相互转换。有阻尼钟摆系统，阻尼项一直耗散能量，系统能量越来越低，最终停止在能量最低的位置。（注：通过增大阻尼项，减少了系统的振荡。在一些机器人系统中，在设计的时候，会改变系统的参，例如阻尼、刚度、质量，来提升系统的性能。）

2.3 引入不变力矩的过阻尼摆

在有阻尼钟摆系统的基础上，如果我们通过电机控制给摆施加一个恒定力矩 \Gamma ，这是我们添加的一个简单的主动控制项，让我们看一下他会怎么影响系统的呢，如下图

系统完整的方程，如下

mL^2\ddot{x}+b\dot{x}+mgLsinx=\mathrm{\Gamma}\\ 式中，

m为质量

L为摆长

b为黏性阻尼系数，阻尼项会阻力 f=b\dot x ，阻碍物体运动

g为重力加速度

\Gamma 为不变力矩（逆时针为正）

在 b 值过大时的过阻尼极限情况下，阻尼项产生的阻力 b\dot x ，要远大于加速度项 mL^2\ddot{x} ，而且在实践过程中，如果系统不是剧烈运动，加速度引入的惯性力相对比较小。为了问题的简化，近似忽略 mL^2\ddot{x} ，则有 b\dot{x}+mgLsinx=\mathrm{\Gamma}\\ 首先，进行无量纲化，两边除以 mgL ，则有 \frac{b}{mgL}\dot{x}=-sinx+\frac{\mathrm{\Gamma}}{mgL}\\ 令 \tau=\frac{mgL}{b}t，r=\frac{\mathrm{\Gamma}}{mgL} ，简化方程 x^\prime=r-sinx, x^\prime=dx/d\tau\\ r为不变力矩和最大重力矩的比值，那么 r ，作用力矩永远不会被平衡掉，单摆将连续的翻转，系统运行轨迹，如下图， r =2,

r =1，不变力矩刚好等于最大重力矩，所有摆会稳定在最高点，夹角180度，运动轨迹，如下图， r=1

r<1 ，不变力矩小于最大重力矩，所以摆会稳定在一个小角度，运动轨迹如下图， r=0.2

注：matlab代码
tspan = [0 150];
y0 = 0;
r = 1;
[t,y] = ode23(@(t,y) (r-sin(y)), tspan, y0);

figure
plot(y*57.3*2,(r-sin(y))*57.3*2)
xlabel('角度 °')
ylabel('角速度 °/s')
title('r=1')

如上三种情况，调节力矩，可以让摆振荡、稳定在最高点或者一个特定的角度。总结来看，通过调节系统参数，例如阻尼，或者，调节控制力矩，可以调节系统运动状态。

控制中主要通过微分方程来对系统建模，研究系统特性的过程就是解微分方程的过程，现代计算机很容易就能够解决求解的方程，例如通过数值计算求解。

几百年前的数学家，把微分方程求解的问题，转化成传递函数，研究特征根（极点）、零点来研究系统特性。通过傅里叶变换，在频域中进行分析，能够通过手动绘制伯德图，就能够快速对系统进行分析。（这些分析工具，都是前人造好的，至于其中的缘由，可能会涉及更高阶的数学理论，没到时机不用强求，先会用就好了。关于其中的理论，建议各位循序渐进的构建自己的架构。）

赠人玫瑰手有余香！记得点赞啦！