1、前言

文章将基于上一篇博客倒立摆自动控制学习笔记（二）的数学模型，设计PID控制器，即

$$\begin{array}{rl}& F-N_{车}-b_1\dot{x}=M\ddot{x}\\ & N_{杆}-f=m(\ddot{x}-\ddot{\theta }lcos\theta +{\dot{\theta }}^{2}lsin\theta )\\ & P_{杆}-mg=-m(\ddot{\theta }lsin\theta +{\dot{\theta }}^{2}lcos\theta )\\ & P_{杆}lsin\theta +N_{杆}lcos\theta -b_2\dot{\theta }=I\ddot{\theta }\\ & N_{车}=N_{杆}\\ & f=\frac{1}{2}{C}_d\rho \dot{x}|\dot{x}|A_{直立}cos\theta \end{array}\text{\hspace{1em}(1-1)}$$

2、级联PID控制器设计

微分方程直接参考上一篇文章倒立摆自动控制学习笔记（三）：LQR控制器设计的结果：
$$\begin{array}{rl}& \ddot{\theta }=19.6\theta -0.4\dot{x}+1.33F\\ & \ddot{x}=3.27\theta -0.27\dot{x}+0.89F\end{array}\text{\hspace{1em}(2-1)}$$
将等号两边拉氏变换，得到
$$\begin{array}{rl}& \Theta (s)=\frac{1.33}{s^2-19.6}F(s)-\frac{0.4s}{s^2-19.6}X(s)\\ & X(s)=\frac{0.89}{s^2+0.27s}F(s)+\frac{3.27}{s^2+0.27s}\Theta (s)\end{array}\text{\hspace{1em}(2-2)}$$
可以得到传递函数
$$\begin{array}{rl}& G_{\theta }(s)=\frac{1.33}{s^2-19.6}\\ & G_x(s)=\frac{0.89}{s^2+0.27s}\end{array}\text{\hspace{1em}(2-3)}$$
注意：这里$\theta$ 和 $x$有很强的耦合，两个传递函数不能看做两个系统。
如果要控制摆角的同时控制小车位移，那么倒立摆就是一个控制输入$F$控制两个控制量 $\theta$ 和 $x$，一般有两种方案：一个外环+一个内环的级联PID控制，以及权重混合的PID控制，在中文互联网能找到的主要是前者，不过后者在谷歌学术找到的论文中采用的更为常见一些。
本文将以级联PID控制方案进行设计，级联PID控制流程图如下：

其中带下标 $d$ 的是给定量。这是摆角控制优先级大于位置控制的情况，也可以反过来（A New Approach on Stabilization Control of an Inverted Pendulum, Using PID Controller）。

内环

我们先看内环。
$-\frac{0.4s}{s^2-19.6}X(s)$视作干扰信号，所以要将其补偿掉。
$$F(s)=PID(s)[{\Theta }_d(s)-\Theta (s)]+0.3sX(s)\text{\hspace{1em}(2-4)}$$

（红色部分为系统固有属性，无法进行任何修改）
显然，系统$G_{\theta }(s)$有一个右半平面极点。
设计PD控制器的$K_p=25$，$K_d=5$，并绘制奈奎斯特图：

% 画出内环（带控制器）开环传递函数的奈奎斯特图、伯德图
clear
clc
f1 = [5 25];%分子
f2 = [0 1.33];%分子
f3 = [0 1];%分母
f4 = [1 0 -19.6];%分母
%conv的用处是把他们求成多项式的形式
num = conv(f1,f2);%求多项式
den = conv(f3,f4);%求多项式
sys = tf(num,den)
%画出所有参数
figure(1)
margin(sys);
% grid on;
%不画，直接给参数赋值出来
%[Gm,Pm,Wcg,Wcp] = margin(sys);%Gm幅值裕度 Pm相角裕度 Wcg穿越频率 Wcp剪切频率
figure(2)
nyquist(sys);  % 画奈奎斯特图
% grid on;

奈奎斯特图绕(-1, j0)点圈数和$G_{\theta }(s)$右半平面极点个数相同，系统稳定。稳定裕度在Bode图中也可以得到。
其实如果想调细一点，可以用相位超前校正来调（根据伯德图计算pid参数方法原理），或者参考《现代控制工程》，直接根据误差系数、相位裕度指标解出一组PID参数，我这里比较偷懒。

外环

这控制框图比我预想中还要复杂一些，第一反应就是上梅森增益公式。但是我们有强大的matlab，用下面的代码可以硬算出蓝色方框里的传递函数。

clear; clc;
syms s x F theta theta_d % 参数

eq1 = F == (theta_d-theta)*(5*s+25) + 0.3*s*x;
eq2 = theta == F*1.33/(s^2-19.6) - x*0.4*s/(s^2-19.6);
eq3 = s*x == theta*3.27/(s+0.27) + F*0.89/(s+0.27);
eqns = [eq1, eq2, eq3];

vars = [x F theta];

% 解方程
sol = solve(eqns, vars);

theta_sol = sol.theta
x_sol = sol.x

% 整理传递函数
G_theta_d_x = x_sol / theta_d;

% 使用collect函数按s的幂次合并同类项
[G_theta_d_x_num, G_theta_d_x_den] = numden(G_theta_d_x);  % 提取分子分母
G_theta_d_x_num_collected = collect(G_theta_d_x_num, s);  % 按s合并分子
G_theta_d_x_den_collected = collect(G_theta_d_x_den, s);  % 按s合并分母
G_theta_d_x = G_theta_d_x_num_collected / G_theta_d_x_den_collected

得到
$$G_{{\theta }_d-x}(s)=\frac{4.45s^3+22.25s^2-65.47s-327.37}{s^4+6.65s^3+13.67s^2+0.02s}\text{\hspace{1em}(2-5)}$$
特征方程
$$D(s)=s(s^3+6.65s^2+13.67s+0.02)\text{\hspace{1em}(2-6)}$$
为简单起见，我们直接看括号里那一部分即可，画劳斯表：

第一列没有出现变号的情况，说明$G_{{\theta }_d-x}(s)$没有右半平面的极点，有一个在虚轴上的极点。
设计控制器
$$K(s)=-\frac{0.01s}{10s+1}\text{\hspace{1em}(2-7)}$$
（控制器看着又怪又四不像，但是我实在搞不出更好的来了）
画出奈奎斯特图：

放大后：

绕(-1, j0)点0圈，系统稳定。
仿真一下看看：

位移有不小的稳态误差，只能说效果很一般。
仿真包括main.m和pendulum_friction_double.slx两个文件，将两个文件放到同一路径，用2024a或更高版本的matlab打开main.m，将代码中的语句：

simout = sim("pendulum_friction_PID", 2);  % 参数为模型名和仿真时间

中的模型名称改为你想要运行的，修改仿真时间，然后运行main.m即可。

仿真文件路径：Pendulum-Control