智能控制——专家控制

⭐⭐⭐✰✰

1 问题描述

考虑一个无人水面艇（USV）航向跟踪控制系统。无人艇的非线性动力学模型如下：
$$
T\ddot{\psi}+KH(\dot{\psi})=K\delta+d(t)
$$
其中：

状态变量$\psi$为航向角(单位：弧度)；$r=\dot{\psi}$为转艏角速度（单位：弧度/秒）
控制输入$\delta$为舵角(单位：弧度)，也是控制器的输出，通常受限，$|\delta|\leq\delta_{max}$
系统参数𝑇为船舶的追随性指数，体现惯性；𝐾为船舶的旋回性指数；$H(\dot{\psi})=\alpha\dot{\psi}+\beta\dot{\psi}^3$ 为一个关键的非线性函数，用于描述船舶在不同转速下的转向特性，参数𝛼、𝛽决定了非线性程度。
外部扰动$d(t)=A_dsin(\omega_dt)+n(t)$，代表综合环境扰动（如风、浪、流），其中 𝑛(t)为随机白噪声。
系统参数：输入舵角$\delta$受限为$|\delta|\leq0.5rad$，$d(t)=0.1sin(0.1t)$模拟持续的环境扰动（如风浪流）。系统参数为$T=2.0s，K=0.8s^{-1}、\alpha=1.0、\beta=0.5$
扰动参数：$A_d=0.1，\omega_d=0.1$，$n(t)$为均值为0，方差为0.01的白噪声
初始条件：$\psi(0)=0.1rad，r(0)=0rad/s$
控制目标：设计控制器，通过调节舵角$\delta$，使船舶航向$\psi$能够跟踪时变的参考航向$\psi_{ref}(t)$（如阶跃信号、斜坡信号或正弦信号）

接下来：

进行专家控制器设计：基于误差 $𝑒=𝜓_{𝑟𝑒𝑓}−𝜓$和误差变化率 $𝑒𝑐=𝑒̇$，设计一个直接映射到舵角𝛿的规则库。
进行专家PID控制器设计：基于(1)中相同的𝑒和𝑒𝑐，设计一个用于动态调整 PID 参数(𝐾𝑝,𝐾𝑖,𝐾𝑑)的规则库，再由标准PID算法计算𝛿
仿真与对比分析：在MATLAB/Simulink中搭建上述船舶模型，在相同条件下仿真实现上述两个控制器，并从超调量、调节时间、稳态误差、控制输入平滑性和鲁棒性（参数摄动下的性能）五个维度进行定量对比

2 问题分析与建模

本题目给出了一个非线性模型，其中控制器输出是δ，被控对象是$ψ$，目的是用专家控制设置控制器和PID控制器来实现$ψ$跟踪$ψ_{ref}$

首先对模型进行分析：
$$
T\ddot{\psi}+KH(\dot{\psi})=K\delta+d(t)
$$
从控制器输出开始，先乘上增益K再加上干扰d(t)得到右边式子，然后需要由此得到$\ddot{\psi}$，再经过两次积分得到$\psi$，因为：
$$
\ddot{\psi}=\frac{1}{T}(K\delta+d(t)-KH(\dot{\psi}))
$$
可以画出：

其中$n(t)$用的Random Number模块作为白噪声；$d(t)$是正弦信号扰动，$\omega_d=0.1，A_d=0.1$；中间的函数模块即$H(\dot{\psi})$；$\psi_0$为初始值0.1rad。前三者具体参数配置如下图。

得到了控制对象的基本模型，就可以把所需要的$e，ec$表示出来然后输入进控制器进行控制了。

3 控制器规则库设计

3.1 专家控制器

专家控制器输入为$e，ec$，输出为δ，用S-Function模块来写，可以得到如图所示的结构。

由S-Function的原理图，我们首先可以写出其主函数：

function [sys,x0,str,ts,simStateCompliance]=expert_control(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);
    otherwise
        DAStudio.error('Simulink:blocks:unhandledFlag', num2str(flag));
end

然后是初始化函数配置相关参数，本控制器没有内部状态变量，输入变量2个，输出变量1个。

function [sys,x0,str,ts,simStateCompliance]=mdlInitializeSizes
sizes = simsizes;
sizes.NumContStates  = 0;%连续状态数
sizes.NumDiscStates  = 0;%离散状态数
sizes.NumOutputs     = 1;%输出量个数
sizes.NumInputs      = 2;%输入量个数
sizes.DirFeedthrough = 1;%是否存在直接馈通。1：存在；0：不存在，缺省为 1 
sizes.NumSampleTimes = 1;%采样时间

sys = simsizes(sizes);
x0  = [];% 设置初始状态
str = [];% 保留变量置空
ts  = [0 0]; % 连续系统
simStateCompliance = 'UnknownSimState';

而对于计算导数、更新状态、计算下一个采样时间、仿真结束四个子函数在此处都用不到，给予默认配置。

而我们的规则库主要体现在计算输出子函数mdlOutputs上

先创建若干条规则：

饱和控制：偏差很大时，打满舵，以求最快速度纠正。
增强控制：偏差中等且仍在变大，用中等舵角抑制其增长。
保守控制：偏差中等但已在缩小，用小舵角柔和纠正，防止超调。
精细控制：接近目标，启用类似PID的精细控制，平滑稳定。
微调控制：默认情况，施加一个微小恒定的纠正力。
输出限制：$|\delta|\leq0.5rad$

再对应转化为代码语言，按照图上的顺序，两个输入中u(1)是e，u(2)是ec

function sys=mdlOutputs(t,x,u)
e=u(1)
ec=u(2);
delta=0;%初始输出为0
if abs(e)>0.5%饱和控制，偏差大时满舵
    delta=0.5*sign(e);
elseif abs(e)>0.2&&ec*e>0%增强控制，偏差中等但是在变大，中等舵角
    delta=0.3*sign(e);
elseif abs(e)>0.2&&ec*e<0%保守控制，偏差中等但是在变小，小舵角
    delta=0.15*sign(e);
elseif abs(e)>0.05%精选控制，接近目标微小舵角
    delta=0.05*e+0.02*ec; 
else
    delta=0.01*e;%默认情况，给一个微小恒定的力
end

delta=min(max(delta,-0.5),0.5);%输出限制
sys=delta;

3.2 专家PID控制器

结构与专家控制器一致，只是控制器的S-Function的规则库不一样。

S-Function大致结构和专家控制器一致，初始化参数设计不同，其规则库体现在状态更新子函数mdlUpdate上

初始化函数：此时定义三个内部状态：累积e积分、上次采样的e、得到的舵角(便于输出)

function [sys,x0,str,ts,simStateCompliance]=mdlInitializeSizes
sizes = simsizes;
sizes.NumContStates  = 3;%连续状态数
sizes.NumDiscStates  = 0;%离散状态数
sizes.NumOutputs     = 1;%输出量个数
sizes.NumInputs      = 2;%输入量个数
sizes.DirFeedthrough = 1;%是否存在直接馈通。1：存在；0：不存在，缺省为 1 
sizes.NumSampleTimes = 1;%采样时间

sys = simsizes(sizes);
x0  = [];% 设置初始状态
str = [];% 保留变量置空
ts  = [0 0]; % 连续系统
simStateCompliance = 'UnknownSimState';

上图为数字PID控制算法，我们所控的Kp，Ki，Kd参数直接就是每个环节前面的系数。

首先设定一组基础PID值，要求能够实现控制，无其他性能指标要求。令$K_P=1.2,K_I=0.3,K_D=0.4$，模拟运行一下：

可见能够实现控制，因此这组默认值是合理的。并由此为基础编写规则库

有若干条规则：

抗饱和：偏差极大，切换为Bang-Bang控制，也就是直接开到满舵，绕过PID。
加速响应：误差在扩大，增强P和D以抑制偏差增长，削弱 I 防积分饱和。
抑制超调：误差在缩小，减弱控制作用，使系统平稳接近目标。
消除静差：接近目标，增强积分作用以消除稳态误差，减弱P和D防振荡。
抗积分饱和：给积分定上下限，同时偏差较大时清零积分。

function sys = mdlUpdate(t,x,u)
Ts=0.01;
e=u(1);
ec=u(2);

integral_e=x(1);%累积积分
pre_e=x(2);%上次的e
delta=x(3);%计算得到的舵角
%基础值
Kp_base=1.2;
Ki_base=0.3;
Kd_base=0.4;
%默认
Kp=Kp_base;
Ki=Ki_base;
Kd=Kd_base;
if abs(e) > 0.5%抗饱和，偏差大直接打满舵，清空积分
    delta=0.5*sign(e);
    integral_e=0;              % 防止积分爆炸
    sys=[integral_e; pre_e; delta];
    return;
end
if e*ec > 0%加速响应，当误差在扩大，增强PD，抑制I
    Kp=Kp_base*1.6;
    Kd=Kd_base*1.6;
    Ki=Ki_base*0.5;
elseif e*ec<0%抑制超调，当误差减小，PID都减小
    Kp=Kp_base*0.7;
    Kd=Kd_base*0.7;
    Ki=Ki_base*0.7;
elseif abs(e)<0.1%消除静差，增强I，抑制PD
    Ki=Ki_base*2.0;
    Kp=Kp_base*0.5;
    Kd=Kd_base*0.5;
end
%更新积分
integral_e=integral_e+e*Ts;
integral_e=max(min(integral_e,1),-1);%抗积分饱和
%求微分
derivative_e=(e-pre_e)/Ts;
%代入PID算法求δ
delta=Kp*e+Ki*integral_e+Kd*derivative_e
delta=max(min(delta,0.5),-0.5);
sys=[integral_e;e;delta];

function sys = mdlOutputs(t,x,u)%输出函数
sys = x(3);

4 模型对比

扰动量初始值等其他参数都相等的情况下，观察控制器控制$ψ$跟踪$ψ_{ref}$d

4.1 阶跃输入

让$ψ_{ref}$为阶跃输入，从0.5阶跃到1

上图左侧为专家控制，右侧为专家PID控制；可以看到当$ψ_{ref}$较小时，控制器受扰动影响较大，输出会产生震荡，并且通过对比，专家PID控制的抗干扰能力比专家控制强，专家控制已经明显受正弦干扰影响而产生周期震荡了。

让$ψ_{ref}$为阶跃输入，从10阶跃到15 。同时为了测试鲁棒性在扰动中新增一个干扰变量。由于δ有限幅，令扰动量为0.2

专家控制响应为：（左侧为阶跃响应ψ，右侧为控制输入δ）

专家PID控制响应为：（左侧为阶跃响应ψ，右侧为控制输入δ）

传统PID控制响应为(参数用专家PID的基础参数$K_P=1.2,K_I=0.3,K_D=0.4$并且设置抗积分饱和)：左边ψ右边δ

由图可得出超调量、调节时间、稳态误差、控制输入的平滑性和鲁棒性

	超调σ	调节时间ts(2%)	稳态误差	δ的平滑性	δ的鲁棒性
专家控制器	2.02%(几乎无超调)	20.882s	0(但是一直在振荡)	δ非常不平滑，各种跳变产生了梳子形状的图象	当给另外的干扰时，δ为适应随之变化，但是最终的响应的均值略微变大，出现了稳态误差，鲁棒性较差
专家PID控制	0.99%(几乎无超调)	15.392s	0(几乎不振荡，较为稳定)	δ虽然也存在锯齿，但是相比于专家控制较为平滑	当给另外的干扰时，δ也随之变化，最终响应不变，没有稳态误差，但是比起扰动前出现了非常轻微的振荡，鲁棒性较好
传统PID控制	1.79%(几乎无超调)	15.196	0(微弱振荡，比专家PID强烈一点，但是整体上还是非常稳定)	δ的平滑性介于专家控制和专家PID控制之间，存在一些特别剧烈变化的锯齿	当给另外的干扰时，δ也随之变化，最终响应不变，没有稳态误差，而且扰动前后的振荡情况几乎完全一致，鲁棒性极强