MATLAB神經網絡（6） PID神經元網絡解耦控制算法——多變量系統控制

6.1 案例背景

6.1.1 PID神經元網絡結構

PID神經元網絡從結構上可以分為輸入層、隱含層和輸出層三層，$n$個控制量的PID神經元網絡包含$n$個並列的相同子網絡，各個子網絡間既相互獨立，又通過網絡連接權值相互聯系。每個子網絡的輸入層有兩個神經元，分別接收控制量的目標值和當前值。每個子網絡的隱含層由比例元、積分元和微分元構成，分別對應着PID控制器中的比例控制、積分控制和微分控制。PID神經元網絡按被控系統控制量的個數可以分為控制單變量系統的單控制量神經元網絡和控制多變量系統的多控制量神經元網絡。其中單控制量神經元網絡是PID神經元網絡的基本形式，多控制量神經元網絡可以看成是多個單控制量神經元網絡的組合形式。

$X_{1}$是控制量的控制目標，$X_{2}$是控制量當前值，$Y$是神經元網絡計算得到的控制律，${\omega _{ij}}$和${\omega _{jk}}$是網絡權值，從中可以看到單控制量神經元網絡是一個三層前向神經元網絡，網絡結構為2-3-1，隱含層包含比例元、積分元和微分元三個神經元。

多控制量神經元網絡可以看成多個單控制量網絡的並聯連接，$X_{11},X_{21},...,X_{n1}$是控制量的控制目標，$X_{12},X_{22},...,X_{n2}$是控制量的當前值，$Y_{1},Y_{2},...,Y_{n}$是多控制量神經元網絡計算得到的控制律，${\omega _{ij}}$和${\omega _{jk}}$是網絡權值。

6.1.2 控制率

• 輸入層：輸出數據$x_{si}$等於輸入數據$X_{si}$：${x_{si}}(k) = {X_{si}}(k)$；
• 隱含層：輸入值\[ne{t_{sj}}(k) = \sum\limits_{i = 1}^2 {{\omega _{ij}}{x_{si}}(k)\;\;\;j = 1,2,3} \]比例神經元\[{u_{s1}} = ne{t_{s1}}(k)\]積分神經元\[{u_{s2}} = ne{t_{s2}}(k) + ne{t_{s2}}(k - 1)\]微分神經元\[{u_{s3}} = ne{t_{s3}}(k) - ne{t_{s3}}(k - 1)\]
• 輸出層：\[{y_k} = \sum\limits_{s = 1}^n {\sum\limits_{j = 1}^3 {{\omega _{jk}}{u_{sj}}(k)} } \]

6.1.3 權值修正

誤差計算公式：\[J = \sum {E = \sum\limits_{k = 1}^n {{{[{y_h}(k) - r(k)]}^2}} } \]其中$n$為輸出節點個數，$y_{h}(k)$為預測輸出，$r(k)$為控制目標。

• 輸出層到隱含層：\[{\omega _{jk}}(k + 1) = {\omega _{jk}}(k) - \eta \frac{{\partial J}}{{\partial {\omega _{jk}}}}\]
• 輸入層到輸出層：\[{\omega _{ij}}(k + 1) = {\omega _{ij}}(k) - \eta \frac{{\partial J}}{{\partial {\omega _{ij}}}}\]

其中$\eta$為學習速率。

6.1.4 控制對象

傳遞函數：

\[\left\{ \begin{array}{l}
{y_1}(k) = 0.4{y_1}(k - 1) + {u_1}(k - 1)/[1 + {u_1}{(k - 1)^2}] + 0.2{u_1}{(k - 1)^3} + 0.5{u_2}(k - 1) + 0.3{y_2}(k - 1)\\
{y_2}(k) = 0.2{y_2}(k - 1) + {u_2}(k - 1)/[1 + {u_2}{(k - 1)^2}] + 0.4{u_2}{(k - 1)^3} + 0.2{u_1}(k - 1) + 0.3{y_3}(k - 1)\\
{y_3}(k) = 0.3{y_3}(k - 1) + {u_3}(k - 1)/[1 + {u_3}{(k - 1)^2}] + 0.4{u_3}{(k - 1)^3} + 0.4{u_2}(k - 1) + 0.3{y_1}(k - 1)
\end{array} \right.\]

6.2 模型建立

$r_{1}$,r_{2},...,r_{n}$是控制器目標，$u_{1}$,u_{2},...,u_{n}$為控制器規律，$y_{1}$,y_{2},...,y_{n}$為控制器當前值。

3個控制量——3個單神經元網絡。

控制量初始值[0,0,0]，控制目標為[0.7,0.4,0.6]，控制時間間隔0.001s。

6.3 編程實現

6.3.1 初始化

%% 基於PID神經網絡的系統控制算法
%% 清空環境變量
clc
clear

%% 網絡結構初始化
rate1=0.006;rate2=0.001; %學習率
K=3;
y_1=zeros(3,1);y_2=y_1;    %輸出值
u_1=zeros(3,1);u_2=u_1;    %控制率
h1i=zeros(3,1);h1i_1=h1i;  %第一個控制量
h2i=zeros(3,1);h2i_1=h2i;  %第二個控制量
h3i=zeros(3,1);h3i_1=h3i;  %第三個控制量
x1i=zeros(3,1);x2i=x1i;x3i=x2i;x1i_1=x1i;x2i_1=x2i;x3i_1=x3i;   %隱含層輸入

%權值初始化
k0=0.03;

%第一層權值
w11=k0*rand(3,2);
w12=k0*rand(3,2);
w13=k0*rand(3,2);
%第二層權值
w21=k0*rand(1,9);
w22=k0*rand(1,9);
w23=k0*rand(1,9);

%值限定
ynmax=1;ynmin=-1;  %系統輸出值限定
xpmax=1;xpmin=-1;  %P節點輸出限定
qimax=1;qimin=-1;  %I節點輸出限定
qdmax=1;qdmin=-1;  %D節點輸出限定
uhmax=1;uhmin=-1;  %輸出結果限定

6.3.2 優化

%% 網絡迭代優化
for k=1:1:300

    %% 控制量輸出計算
    %--------------------------------網絡前向計算--------------------------
    
    %系統輸出
    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.4*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(yn>ynmax)=ynmax;
    yn(yn<ynmin)=ynmin;
    
    %輸入層輸出
    x1o=[r1(k);yn(1)];x2o=[r2(k);yn(2)];x3o=[r3(k);yn(3)];
    
    %隱含層 
    x1i=w11*x1o;
    x2i=w12*x2o;
    x3i=w13*x3o;

    %比例神經元P計算
    xp=[x1i(1),x2i(1),x3i(1)];
    xp(xp>xpmax)=xpmax;
    xp(xp<xpmin)=xpmin;
    qp=xp;
    h1i(1)=qp(1);h2i(1)=qp(2);h3i(1)=qp(3);

    %積分神經元I計算
    xi=[x1i(2),x2i(2),x3i(2)];
    qi_1=[h1i(2),h2i(2),h3i(2)];
    qi=qi_1+xi;
    qi(qi>qimax)=qimax;
    qi(qi<qimin)=qimin;
    h1i(2)=qi(1);h2i(2)=qi(2);h3i(2)=qi(3);

    %微分神經元D計算
    xd=[x1i(3),x2i(3),x3i(3)];
    xd_1=[x1i_1(3),x2i_1(3),x3i_1(3)];
    qd=xd-xd_1;
    qd(qd>qdmax)=qdmax;
    qd(qd<qdmin)=qdmin;
    h1i(3)=qd(1);h2i(3)=qd(2);h3i(3)=qd(3);

    %輸出層計算
    wo=[w21;w22;w23];
    qo=[h1i',h2i',h3i'];qo=qo';
    uh=wo*qo;
    uh(uh>uhmax)=uhmax;
    uh(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)=(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;
    
    %調整w22
    Sig2=sign(ypc./(uh(2)+0.00001));
    dw22=sum(error.*Sig2)*qo';
    w22=w22+rate2*dw22;
    
    %調整w23
    Sig3=sign(ypc./(uh(3)+0.00001));
    dw23=sum(error.*Sig3)*qo';
    w23=w23+rate2*dw23;

    %輸入層和隱含層權值調整
    delta2=zeros(3,3);
    wshi=[w21;w22;w23];
    %調整w11
    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.00001));
    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
        dws1(p,:)=alter(p)*s1(p,:);
    end
    w11=w11+rate1*dws1;

    %調整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
            alter2(j)=alter2(j)+delta2(p,:)*wshi2_2(:,j);
        end
    end
    for p=1:1:3
        dws2(p,:)=alter2(p)*s2(p,:);
    end
    w12=w12+rate1*dws2;
    
    %調整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;

    %參數更新
    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;
end

6.3.3 結果分析

%% 結果分析
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);
legend('控制目標','實際輸出');

subplot(3,1,2)
plot(time,r2,'r-',time,y2,'b-');

ylabel('控制量2','fontsize',12);
legend('控制目標','實際輸出');
subplot(3,1,3)
plot(time,r3,'r-',time,y3,'b-');       
xlabel('時間（秒）','fontsize',12);ylabel('控制量3','fontsize',12);
legend('控制目標','實際輸出');

figure(2)
plot(time,u1,'r-',time,u2,'g-',time,u3,'b');
title('PID神經網絡提供給對象的控制輸入');
xlabel('時間'),ylabel('被控量');
legend('u1','u2','u3');grid

figure(3)
figure(3)
plot(time,J,'r-');
axis([0,0.2,0,1]);grid
title('控制誤差曲線','fontsize',12);
xlabel('時間','fontsize',12);ylabel('控制誤差','fontsize',12);

6.3.4 Simulink仿真

6.4 擴展

6.4.1 增加動量項

\[\begin{array}{l}
{\omega _{jk}}(k + 1) = {\omega _{jk}}(k) - \eta \frac{{\partial J}}{{\partial {\omega _{jk}}}} + {\eta _1}[{\omega _{jk}}(k) - {\omega _{jk}}(k - 1)]\\
{\omega _{ij}}(k + 1) = {\omega _{ij}}(k) - \eta \frac{{\partial J}}{{\partial {\omega _{ij}}}} + {\eta _1}[{\omega _{ij}}(k) - {\omega _{ij}}(k - 1)]
\end{array}\]

6.4.2 神經元系數

\[\begin{array}{l}
{u_{s1}} = {k_p} \cdot ne{t_{s1}}(k)\\
{u_{s2}} = {k_i} \cdot [ne{t_{s2}}(k) + ne{t_{s2}}(k - 1)]\\
{u_{s3}} = {k_d} \cdot [ne{t_{s3}}(k) - ne{t_{s3}}(k - 1)]
\end{array}\]

6.4.3 PID神經元網絡權值優化

粒子群/遺傳算法——解決可能陷入局部最優的問題，改善控制效果。