PIDNN

博主主要做汽车控制方向，之前也学了一些Neural Network相关的东西，在Tensorflow里面还是蛮好搭一个小的NN来玩一下的。印象里面这些神经网格算法总是和感知、认知等方向关联，最近正好在看清华大学Matlab智能算法一书，里面的PIDNN算法给了我很大的启发，其自适应性能帮助工程师减少很多繁杂的标定工作

原理上看如下图所示，本质上还是输入的前向传播和误差反向传播。

把书里面的代码整理一下，便于后续修改使用，我后面也会慢慢增加相关代码注释，慢慢看其实并不复杂，需要耐下心来慢慢学习。代码运行在Matlab里，没有使用复杂的库函数，正常都能运行。

clc
clear all

% 网格初始化
rate1 = 0.001;
rate2 = 0.005;
rate3 = 0.0001;
k = 0.1;
K = 5;
y_1 = zeros(3,1);
y_2 = y_1;
y_3 = y_2;
u_1 = zeros(3,1);
u_2 = u_1;
u_3 = u_2;  % 控制率

h1i = zeros(3,1);
h1i_1 = h1i; % 第一个控制量   Dim 3*1
h2i = zeros(3,1);
h2i_1 = h2i; % 第二个控制量
h3i = zeros(3,1);
h3i_1 = h3i; % 第三个控制量
x1i = zeros(3,1);
x2i = x1i;
x3i = x2i;
x3i = x2i;
x1i_1 = x1i;
x2i_1 = x2i;
x3i_1 = x3i; % 隐含层输出

% 权值初始化
k0 = 0.03;

% 第一层权值
w11 = k0 * rand(3,2);w11_1 = w11; w11_2 = w11_1; % Dim 3*2
w12 = k0 * rand(3,2);w12_1 = w12; w12_2 = w12_1;
w13 = k0 * rand(3,2);w13_1 = w13; w13_2 = w13_1;

% 第二层权值
w21 = k0 * rand(1,9);w21_1 = w21; w21_2 = w21_1;  % Dim 1*9
w22 = k0 * rand(1,9);w22_1 = w22; w22_2 = w22_1;
w23 = k0 * rand(1,9);w23_1 = w23; w23_2 = w23_1;

% 值限定
ynmax = 1; ynmin = -1; % 系统输出值限定
xpmax = 1; xpmin = -1; % P节点输出限定
qimax = 1; qimin = -1; % I节点输出限定
qdmax = 1; qdmin = -1; % D节点输出限定
uhmax = 1; uhmin = -1; % 输出结果限定
%% 网络迭代优化
for k = 1:1:200
    % 控制量输出计算
    % 网格前向计算
    
    % 系统输出
    y1(k) = (0.4 * y_1(1) + u_1(1)/(1+u_1(1)^2 )+0.2*u_1(1)^3+0.5*u_1(2))+0.3*y_1(2);
    y2(k) = (0.2 * y_1(2) + u_1(2)/(1+u_1(2)^2) +0.4*u_1(2)^3+0.2*u_1(1))+0.3*y_1(3);
    y3(k) = (0.3 * y_1(3) + u_1(3)/(1+u_1(3)^2 )+0.4*u_1(3)^3+0.2*u_1(2))+0.3*y_1(1);
    
    % 控制目标
    r1(k) = 0.7;r2(k) = 0.4; r3(k) =0.6;
    
    % 系统输出限制
    yn = [y1(k), y2(k), y3(k)];
    yn(find(yn>ynmax)) = ynmax;
    yn(find(yn<ynmin)) = ynmin;
    
    % 输入层输出
    x1o = [r1(k);yn(1)]; % Dim 2*1
    x2o = [r2(k);yn(2)];
    x3o = [r3(k);yn(3)];
    
    % 隐含层
    x1i = w11*x1o; % Dim3*2  * Dim 2*1 = Dimi 3*1
    x2i = w12*x2o;
    x3i = w13*x3o;
    
    % 比例神经元P计算
    xp = [x1i(1),x2i(1),x3i(1)]; % Dim 1*3 由于有三个输入，对应每个值都有其比例系数，积分系数和微分系数
    xp(find(xp > xpmax)) = xpmax;
    xp(find(xp > xpmax)) = xpmax;
    qp = xp;
    h1i(1) =qp(1);  h2i(1) =qp(2); h3i(1) =qp(3);
    
    % 积分神经元I计算
    xi = [x1i(2),x2i(2),x3i(2)]; % Dim 1*3
    qi = [0,0,0];
    qi_1 = [h1i(2),h2i(2),h3i(2)];% 上一个循环的值
    qi = qi_1 + xi; % 上一个循环的值+本循环的值 = 积分
    qi(find(qi>qimax)) = qimax; % 限制最大
    qi(find(qi<qimin)) = qimin; % 限制最小
    h1i(2) = qi(1); h2i(2) = qi(2); h3i(2) =qi(3);
    
    % 微分神经元D计算
    xd = [x1i(3),x2i(3),x3i(3)]; % Dim 1*3
    qd = [0,0,0];
    xd_1 = [x1i_1(3),x2i_1(3),x3i_1(3)]; % 上一个循环的值
    qd = xd - xd_1; % 本循环-上个循环的值 = 微分
    qd(find(qd>qimax)) = qimax;
    qd(find(qd<qimin)) = qimin;
    h1i(3) = qd(1); h2i(3) = qd(2); h3i(3) =qd(3);
    
    % 输出层计算
    wo = [w21;w22;w23]; % Dim 3*9
    qo = [h1i',h2i',h3i']; % Dim 1*9 
    qo = qo';
    uh = wo * qo; % Dim3*1 输出结果
    uh(find(uh>uhmax)) = uhmax;
    uh(find(uh<uhmin)) = uhmin;
    u1(k) = uh(1);
    u2(k) = uh(2);
    u3(k) = uh(3);
    
    % 网络反馈修正
    % 计算误差
    error = [r1(k) - y1(k);r2(k) - y2(k);r3(k) - y3(k)]; % 求系统与真实误差
    error1(k) = error(1);error2(k) = error(2); error3(k) = error(3);
    J(k) =0.5*(error(1)^2+ error(2)^2 +error(3)^2); % 构造损失函数
    ypc = [y1(k)-y_1(1);y2(k) - y_1(2);y3(k)-y_1(3)];
    uhc = [u_1(1)-u_2(1);u_1(2)-u_2(2);u_1(3)-u_2(3)];
    
    % 隐含层和输出层权值调整
    
    % 调整 w21
    Sig1 = sign(ypc./(uhc(1) + 0.00001));
    dw21 = sum(error .* Sig1)*qo';
    w21 = w21+rate2*dw21 +rate3*(w21_1 - w21_2);
    
    % 调整w22
    Sig2 = sign(ypc./(uh(2) + 0.00001));
    dw22 = sum(error .* Sig2)*qo';
    w22 = w22+rate2*dw22 +rate3*(w22_1 - w21_2);
    
    % 调整w23
    Sig3 = sign(ypc./(uh(3) + 0.00001));
    dw23 = sum(error .* Sig3)*qo';
    w23 = w23+rate2*dw23 +rate3*(w23_1 - w23_2);
    
    % 输入层和隐含层权值调整
    delta2 = zeros(3,3);
    wshi = [w21;w22;w23];
    for t = 1:1:3
        delta2(1:3,t) = error(1:3).*sign(ypc(1:3)./(uhc(t)+0.00000001));
    end
    for j = 1:1:3
        sgn(j) = sign((h1i(j)-h1i_1(j))/(x1i(j) -x1i_1(j)+0.0000001));
    end
    
    s1 = sgn'*[r1(k),y1(k)];
    wshi2_1 = wshi(1:3,1:3);
    alter =zeros(3,1);
    dws1 = zeros(3,2);
    for j = 1:1:3
        for p = 1:1:3
            alter(j) = alter(j) + delta2(p,:)*wshi2_1(:,j);
        end
    end
    
    for p = 1:1:3
        dwsl(p,:) = alter(p)*s1(p,:);
    end
    w11 = w11 + rate1 *dws1 + rate3*(w11_1 - w11_2);
    
    % 调整w12
    for j =1:1:3
        sgn(j) = sign((h2i(j)-h2i_1(j))/(x2i(j) -x2i_1(j)+0.0000001));
    end
    s2 =  sgn'*[r2(k),y2(k)];
    wshi2_2 = wshi(:,4:6);
    alter2 =zeros(3,1);
    dws2 = zeros(3,2);
    for j =1:1:3
        for p = 1:1:3
             alter(j) = alter(j) + delta2(p,:)*wshi2_2(:,j);
        end
    end
    for p = 1:1:3
        dws2(p,:) = alter2(p)*s2(p,:);
    end
    w12 = w12 + rate1 *dws2 + rate3*(w12_1 - w12_2);
    
    % 调整w13
     for j =1:1:3
        sgn(j) = sign((h3i(j)-h3i_1(j))/(x3i(j) -x3i_1(j)+0.0000001));
    end
    s3 =  sgn'*[r3(k),y3(k)];
    wshi2_3 = wshi(:,7:9);
    alter3 =zeros(3,1);
    dws3 = zeros(3,2);
    for j =1:1:3
        for p = 1:1:3
             alter3(j) = alter3(j) + delta2(p,:)*wshi2_3(:,j);
        end
    end
    for p = 1:1:3
        dws3(p,:) = alter2(p)*s3(p,:);
    end
    w13 = w13 + rate1 *dws3 + rate3*(w13_1 - w13_2);
    
    % 参数更新
    u_3 = u_2; u_2 = u_1; u_1 = uh;
    y_2 =y_1; y_1 = yn;
    h1i_1 = h1i; h2i_1 = h2i; h3i_1 = h3i;
    x1i_1 = x1i; x2i_1 = x2i; x3i_1 = x3i;
    w11_1 = w11; w11_2 = w11_1;
    w12_1 = w12; w12_2 = w12_1;
    w13_1 = w13; w13_2 = w13_1;
    % 第二层权值
    w21_1 = w21; w21_2 = w21_1;
    w22_1 = w22; w22_2 = w22_1;
    w23_1 = w23; w23_2 = w23_1;
end

    % 结果分析
    time = 0.001 *(1:k);
    figure(1)
    subplot(3,1,1)
    plot(time,r1,'r-',time,y1,'b-');
    title('PID神经元网络控制','fontsize',12);
    ylabel('控制量1','fontsize',12);
    
    subplot(3,1,2)
    plot(time,r2,'r-',time,y2,'b-');
    ylabel('控制量2','fontsize',12);
    
    subplot(3,1,3)
    plot(time,r3,'r-',time,y3,'b-');
    xlabel('Time(s)','fontsize',12);
    ylabel('控制量3','fontsize',12);

    figure(2)
    plot(time,u1,'r-',time,u2,'g-',time,u3,'b');
    titile('PID神经网络提供给对象的控制输入');
    xlabel('时间'),ylabel('被控量');
    legend('u1','u2','u3');grid
    figure(3)
    plot(time,J,'r-');
    axis([0 ,0.2,0.1]);grid
    title('被控误差曲线','fontsize',12);
    xlabel('时间','fontsize',12);
    xlabel('时间','fontsize',12);
    ylabel('被控误差','fontsize',12);