智能控制——神经网络控制

1 问题描述

考虑与前面作业中相同的无人水面艇（USV）航向跟踪控制系统，通过调节舵角$δ$，使船舶航向$ψ$能够跟踪时变的参考航向$𝜓_{𝑟𝑒𝑓}(t)$，参考信号为：正弦信号$\psi_{ref}(t)=0.2sin(0.1t)$

(1)在matlab中利用已给的带扰动噪声系统模型采集系统的输入输出数据，采用BP神经网络作为辨识模型，建立USV的神经网络辨识模型，并分析辨识结果的精度。
(2)设计一个BP神经网络直接作为控制器，画出该直接神经网络控制系统的结构框图，利用梯度下降法，推导该神经网络控制器的权值在线更新律，并进行控制系统仿真，在 仿真中分别符号近似法和神经网络估计法估计未知系统Jacobian。
(3)用一个BP神经网络在线调整PID控制器的参数，画出该神经网络PID控制系统的结构框图，同样利用梯度下降法，推导神经网络的权值更新律，并进行控制系统仿真。
(4)借鉴模糊自适应相同的方式，设计针对仿射非线性USV模型的神经网络自适应控制器，给出控制律和自适应律，并进行仿真。（提示，模糊逼近可以用RBF神经网络逼近代替，RBF网络即为激发函数为高斯函数的前向网络）

2 神经网络辨识模型

$$
T\ddot{\psi}+KH(\dot{\psi})=K\delta+d(t)
$$
代入参数
$$
2\ddot{\psi}+0.8(\dot{\psi}+0.5\dot{\psi}^3)=0.8\delta+d(t)\
$$
为方便训练将模型搭建为如下形式

选用NN Predict Controller模块的Plant Identification

参数设置如图所示，训练数据输入δ范围限制为[-0.5，0.5]，输出暂限制为[-10，10]

随机输入信号收集8000个样本

训练结果如图

辨识精度较高

预测控制结果为：

3 神经网络控制器

3.1 控制框图和控制律推导

控制器系统框图为

此处用的网络是三层3-5-1前馈BP网络：输入层$N_i=3$，隐藏层$N_h=5$，输出层$N_o=1$，激活函数选tanh

输出层梯度推导:用jacobian有$ \frac{\partial y}{\partial u}=J$，等效输出层误差信号$\delta_o(k)=-\frac{e(k)}{J(k)}$

性能指标：
$$
E(k) = \frac{1}{2} e^2(k)
$$
梯度下降：
$$
\Delta \mathbf{W}_2(k) = -\eta \frac{\partial E}{\partial \mathbf{W}_2}
$$
链式法则：
$$
\frac{\partial E}{\partial \mathbf{W}_2} = \frac{\partial E}{\partial u} \frac{\partial u}{\partial \mathbf{W}_2}
$$
其中：
$$
\frac{\partial E}{\partial u} = \frac{\partial E}{\partial e} \frac{\partial e}{\partial y} \frac{\partial y}{\partial u} = e(k)(-1)J(k)\
\frac{\partial u}{\partial W_2}=h^T(k)
$$
输出层权值更新律：
$$
\Delta \mathbf{W}_2(k) = \eta,\delta_o(k),\mathbf{h}^T(k) + \alpha,\Delta \mathbf{W}_2(k-1) \ \mathbf{W}_2(k+1) = \mathbf{W}_2(k) + \Delta \mathbf{W}_2(k)
$$
隐藏层误差信号为：
$$
\delta_h(k) = \left(\mathbf{W}_2^T(k)\delta_o(k)\right) \odot \left(1 - \mathbf{h}^2(k)\right)
$$
隐藏层权值更新律:
$$
\Delta \mathbf{W}_1(k) = \eta,\delta_h(k),\mathbf{x}^T(k) + \alpha,\Delta \mathbf{W}_1(k-1) \\ \mathbf{W}_1(k+1) = \mathbf{W}_1(k) + \Delta \mathbf{W}_1(k)
$$

3.2 仿真建模

由结果框图搭建出如下仿真模型

其中USV就是整个被控模型：

其中BP控制器代码为

function u = BP_Controller(ek, ek1, uk1, J, step)
eta       = 0.02;
alpha     = 0.7; 
u_max     = 0.5; 
u_min     = -0.5;
Ni = 3;
Nh = 5; 
No = 1; 
persistent W1 W2 dW1 dW2
if isempty(W1)
    W1  = randn(Nh,Ni)*0.2;
    W2  = randn(No,Nh)*0.2;
    dW1 = zeros(size(W1));
    dW2 = zeros(size(W2));
end
x = [ek; ek1; uk1];
net_h = W1*x;
h     = tanh(net_h);
u     = W2*h;
if step < 5
    u = u + 0.01*randn(); %先小随机输出
end
u = max(min(u, u_max), u_min);
if step > 1
    J = max(abs(J), 1e-3) * sign(J);
    % 反向传播
    delta_o = -ek / J;
    %输出层更新
    dW2_new = eta * delta_o * h' + alpha * dW2; 
    %隐藏层更新
    delta_h = (W2' * delta_o) .* (1 - h.^2);
    dW1_new = eta * delta_h * x' + alpha * dW1;
    %权值更新
    W2 = W2 + dW2_new;
    W1 = W1 + dW1_new;
    dW2 = dW2_new;
    dW1 = dW1_new;
end
end

jacobian的计算方法可有符号近似法和神经网络估计法

3.3 运行结果

3.3.1 符号近似法

function J = Jacobian(yk, yk1, uk, uk1)
dy = yk - yk1;
du = uk - uk1;
if abs(du) < 1e-3
    du = sign(du) * 1e-3 + 1e-3;
end
J = dy / du;
end

(左图为纯输出，右图为输出与参考信号对比)

3.3.2 神经网络近似法

function Jhat = NN_Jacobian(yk, yk1, uk, uk1)
eta = 0.02;  
Ni = 2;  %输入层
Nh = 4;  %隐藏层
No = 1;  %输出层
persistent W1 W2 dW1 dW2
if isempty(W1)
    W1 = randn(Nh,Ni)*0.1;  
    W2 = randn(No,Nh)*0.1; 
    dW1 = zeros(size(W1)); 
    dW2 = zeros(size(W2));
end
dy = yk - yk1;
du = uk - uk1;
if abs(du) < 1e-3
    du = sign(du)*1e-3 + 1e-3; %防止分母为0
end
J_real = dy / du;
%前向计算
x = [yk; uk];
h = tanh(W1*x);
Jhat = W2 * h;
%反向传播
eJ = J_real - Jhat;
delta_o = eJ;
delta_h = (W2' * delta_o) .* (1 - h.^2);
%权值更新
dW2_new = eta * delta_o * h' + 0.8*dW2;
dW1_new = eta * delta_h * x' + 0.8*dW1;
W2 = W2 + dW2_new;
W1 = W1 + dW1_new;
dW2 = dW2_new;
dW1 = dW1_new;
end

(左图为纯输出，右图为输出与参考信号对比)

用两种方法计算jacobian结果是一样的

4 神经网络PID控制器

4.1 结果框图和控制律

此处用的网络是三层3-5-1前馈BP网络：输入层$N_i=3$，隐藏层$N_h=6$，输出层$N_o=3$，激活函数选tanh

输出层梯度推导：令$J=e^2+0.05u^2$，$\delta_o= - \begin{bmatrix} J\ J\ J \end{bmatrix}$
性能指标：
$$
J(k) = e^2(k) + \lambda u^2(k), \quad \lambda = 0.05
$$
梯度下降：
$$
\mathbf{W}_2 \leftarrow \mathbf{W}_2 - \eta \frac{\partial J}{\partial \mathbf{W}_2}
$$
链式法则
$$
\frac{\partial J}{\partial \mathbf{W}_2} = \frac{\partial J}{\partial \Delta \mathbf{K}} \frac{\partial \Delta \mathbf{K}}{\partial \mathbf{W}_2}
$$
其中：
$$
\Delta \mathbf{K} = \mathbf{W}_2 \mathbf{h} \Rightarrow \frac{\partial \Delta \mathbf{K}}{\partial \mathbf{W}_2} = \mathbf{h}^T
$$
输出层权值更新律：
$$
\Delta \mathbf{W}_2(k) = \eta,\delta_o(k),\mathbf{h}^T(k) + \alpha,\Delta \mathbf{W}_2(k-1) \ \mathbf{W}_2(k+1) = \mathbf{W}_2(k) + \Delta \mathbf{W}_2(k)
$$
隐藏层误差信号为：
$$
\delta_h(k) = \left(\mathbf{W}_2^T(k)\delta_o(k)\right) \odot \left(1 - \mathbf{h}^2(k)\right)
$$
隐藏层权值更新律:
$$
\Delta \mathbf{W}_1(k) = \eta,\delta_h(k),\mathbf{x}^T(k) + \alpha,\Delta \mathbf{W}_1(k-1) \\ \mathbf{W}_1(k+1) = \mathbf{W}_1(k) + \Delta \mathbf{W}_1(k)
$$

4.2 仿真建模及结果

由系统框图得到如下建模

其中Ts是通过clock模块加unit delay模块计算得到的，step是通过counter limited模块计算得到的
PID模块代码为：

function [u, u_int_next] = PID_A(Kp, Ki, Kd, e, e1, u_int_prev, Ts)
de = (e - e1)/Ts;
if abs(u_int_prev) < 1
    u_int = u_int_prev + e*Ts;
else
    u_int = u_int_prev;
end
%控制律
u = Kp*e + Ki*u_int + Kd*de;
%限幅
u = max(min(u,0.5), -0.5);
u_int_next = u;
end

BP_PID整定神经网络代码为：

function [Kp_out, Ki_out, Kd_out] = BP_Controller_PID_A(e, e1, u, step)
Nh = 6;%隐藏层
Ni = 3;%输入层
No = 3;%输出层
eta = 0.002;
alpha = 0.6;
N_update = 20;
%初始PID值
persistent Kp Ki Kd
if isempty(Kp)
    Kp = 1.2;
    Ki = 0.3;
    Kd = 0.05;
end
persistent W1 W2 dW1 dW2
if isempty(W1)
    W1  = 0.3 * randn(Nh,Ni);
    W2  = 0.2 * randn(No,Nh);
    dW1 = zeros(size(W1));
    dW2 = zeros(size(W2));
en
de = e - e1;
x  = [e; de; u];
%前向计算
h  = tanh(W1 * x);
dK = W2 * h;
dKp = max(min(dK(1), 0.01), -0.01);
dKi = max(min(dK(2), 0.005), -0.005);
dKd = max(min(dK(3), 0.002), -0.002);

Kp = max(Kp + dKp, 0);
Ki = max(Ki + dKi, 0);
Kd = max(Kd + dKd, 0);

Kp = min(Kp, 5);
Ki = min(Ki, 2);
Kd = min(Kd, 0.2);

sigma = 1e-4;
Kp = Kp - sigma * (Kp - 1.2);
Ki = Ki - sigma * (Ki - 0.3);
Kd = Kd - sigma * (Kd - 0.05);

if mod(step, N_update) == 0
    J = e^2 + 0.05*u^2;

    delta_o = -J * ones(No,1);

    dW2_new = eta * delta_o * h' + alpha * dW2;
    delta_h = (W2' * delta_o) .* (1 - h.^2);
    dW1_new = eta * delta_h * x' + alpha * dW1;

    W2 = W2 + dW2_new;
    W1 = W1 + dW1_new;

    dW2 = dW2_new;
    dW1 = dW1_new;
end
Kp_out = Kp;
Ki_out = Ki;
Kd_out = Kd;
end

得到的结果为：(左图为纯输出，右图为输出与参考信号对比)

效果比神经网络直接控制更好，但是仍然存在大量波动。

5 神经网络自适应控制器

5.1 控制律和自适应律

仿射非线性形式：令$x_1=\psi$，$\dot{x_1}=x_2=\dot{\psi}$，$\dot{x_2}=0.4x_2+0.2x_2^3+0.4u$
$$
\ddot{\psi}=0.4\dot{\psi}+0.2\dot{\psi}^3+0.4\delta\
\begin{cases} \dot{x}_1 = x_2 \ \dot{x}_2 = 0.4x_2 + 0.2x_2^3 + 0.4u \end{cases}
$$
结果框图为：

假设系统为未知非线性系统$\ddot{\psi} = f(x_2) + g(x_2),u$，使用RBF网络来逼近$\hat f(x_2) \approx f(x_2),\quad \hat g(x_2) \approx g(x_2)$
RBF神经网络结构为：
$$
\Phi_i(x_2) = \exp!\left( -\frac{(x_2 - c_i)^2}{\sigma^2} \right)\\Phi(x_2) = [\Phi_1,\dots,\Phi_N]^T\
\begin{aligned} \hat f(x_2) &= W_f^T \Phi(x_2) \ \hat g(x_2) &= W_g^T \Phi(x_2),\quad \hat g \ge g_{\min} > 0 \end{aligned}
$$
定义误差与滑模变量$e=\psi-\psi_d，s=\dot{e}+\lambda e$

控制律：
$$
u = \frac{ -\hat f(x_2) - k,s + \ddot{\psi}_d - \lambda e_2 }{ \hat g(x_2) }\
|u| \le 0.5
$$
自适应律：
$$
\dot W_f = -\gamma_f, s, \Phi(x_2) - \sigma_f W_f \W_f(k+1) = W_f(k) + \dot W_f\
\dot W_g = \gamma_g, s, u, \Phi(x_2) - \sigma_g W_g
$$

5.2 仿真建模和结果

由此建模

其中给Wf设置初始值为zeros(5,1)；Wg设置初始值为ones(5,1)*0.5

误差处理模块，用于计算$e,\dot{e},s$

function [e2, s] = Error_Block(x1, x2, psi_d, psi_d_dot)
lambda = 1.5;
e1 = x1 - psi_d;
e2 = x2 - psi_d_dot;
s = e2 + lambda * e1;
end

RBF逼近

function [f_hat, g_hat, Phi] = RBF_NN(x2, Wf, Wg, f_hat_prev, g_hat_prev, s)
N = length(Wf);
c = linspace(-2,2,N);
sigma = 1.0;
x2_clip = min(max(x2, c(1)), c(end));
%基函数
Phi = zeros(N,1,'double');
for i = 1:N
    Phi(i) = exp(-(x2_clip - c(i))^2 / sigma^2);
end
%NN输出
f_hat_raw = sum(Wf(:) .* Phi(:));
g_hat_raw = max(sum(Wg(:) .* Phi(:)), 0.1);
%动态低通滤波
s_threshold = 0.5;
if abs(s) > s_threshold
    alpha = 0.9; 
else
    alpha = 1.0; 
end
f_hat = alpha * f_hat_prev + (1-alpha) * f_hat_raw;
g_hat = alpha * g_hat_prev + (1-alpha) * g_hat_raw;
end

控制律

function u = Control_Law(f_hat, g_hat, s, e2, psi_d_ddot)
lambda = 1.5;
k = 2.0;
g_hat = max(g_hat, 0.05);
u_raw = (-f_hat - k*s + psi_d_ddot - lambda*e2) / g_hat;
Umax = 0.5;
u = min(max(u_raw, -Umax), Umax);
end

自适应律

function [Wf_next, Wg_next] = Adaptive_Law(s, u, Phi, Wf, Wg)
gamma_f = 0.5;
gamma_g = 0.2;
%更新律
Wf_dot = -gamma_f * s * Phi - 0.001*Wf;
Wg_dot =  gamma_g * s * u * Phi - 0.001*Wg;
%更新权值
Wf_next = Wf + Wf_dot;
Wg_next = Wg + Wg_dot;
%限幅
Wf_max = 10;
Wg_max = 10;
Wg_min = 0.1;
Wf_next = min(max(Wf_next, -Wf_max), Wf_max);
Wg_next = min(max(Wg_next, Wg_min), Wg_max);
end

被控对象USV模型

function [x1_dot,x2_dot]= USV(u,x1,x2)
x1_dot = x2;
x2_dot = 0.4*x2 + 0.2*x2^3 + 0.4*u;
end

得到的结果为：(左图为纯输出，右图为输出与参考信号对比)

神经网络自适应控制器的效果非常好，几乎0误差而且时间长也不会出现失调现象。

在控制器输出的u上加白噪声扰动再运行，发现几乎依旧是0误差，鲁棒性非常强。