基于神经网络的预测控制模型仿真.doc

基于神经网络的预测控制模型仿真摘要：本文利用一种权值可以在线调整的动态BP神经网络对模型预测误差进行拟合并与预测模型一起构成动态组合预测器，在此基础上形成对模型误差具有动态补偿能力的预测控制算法。该算法显著提高了预测精度，增强了预测控制算法的鲁棒性。关键词：预测控制 神经网络 动态矩阵 误差补偿1.引言动态矩阵控制(DMC)是一种适用于渐近稳定的线性或弱非线性对象的预测控制算法， 目前已广泛应用于工业过程控制。它基于对象阶跃响应系数建立预测模型， 因此建模简单， 同时采用多步滚动优化与反馈校正相结合，能直接处理大时滞对象，并具有良好的跟踪性能和较强的鲁棒性。 但是，DMC算法在实际控制中存在一系列问题，模型失配是其中普遍存在的一个问题，并会不同程度地影响系统性能。DMC在实际控制中产生模型失配的原因主要有2个，一是诸如建模误差、环境干扰等因素，它会在实际控制的全程范围内引起DMC的模型失配；二是实际系统的非线性特性，这一特性使得被控对象的模型发生变化，此时若用一组固定的阶跃响应数据设计控制器进行全程范围的控制，必然会使实际控制在对象的非建模区段内出现模型失配。针对DMC模型失配问题，已有学者进行了大量的研究，并取得了丰富的研究成果，其中有基于DMC控制参数在线辨识的智能控制算法，基于模型在线辨识的自校正控制算法以及用神经元网络进行模型辨识、在辨识的基础上再进行动态矩阵控制等。这些算法尽管进行在线辨识修正对象模型参数，仍对对象降阶建模误差(结构性建模误差)的鲁棒性不好，并对随机噪声干扰较敏感。针对以上问题，出现了基于误差校正的动态矩阵控制算法。这些文献用基于时间序列预测的数学模型误差代替原模型误差，得到对未来误差的预测。有人还将这种误差预测方法引入动态矩阵控制，并应用于实际。这种方法虽然使系统表现出良好的稳定性，但建立精确的误差数学模型还存在一定的困难。 本文利用神经网络通过训练学习能逼近任意连续有界函数的特点，建立了一种采用BP神经网络进行预测误差补偿的DMC预测控制模型。其中神经网络预测误差描述了在预测模型中未能包含的一切不确定性信息，可以归结为用BP神经网络基于一系列过去的误差信息预测未来的误差，它作为模型预测的重要补充，不仅降低建立数学模型的负担，而且还可以弥补在对象模型中已简化或无法加以考虑的一切其他因素。 本文通过进行仿真，验证了基于神经网络误差补偿的预测控制算法的有效性及优越性，从而增强了动态矩阵控制算法的鲁棒性。2.BP神经网络模型BP神经网络模型处理信息的基本原理是：输入信号Xi通过中间节点（隐层点）作用于输出节点，经过非线形变换，产生输出信号Yk，网络训练的每个样本包括输入向量X和期望输出量t，网络输出值Y与期望输出值t之间的偏差，通过调整输入节点与隐层节点的联接强度取值Wij和隐层节点与输出节点之间的联接强度Tjk以及阈值，使误差沿梯度方向下降，经过反复学习训练，确定与最小误差相对应的网络参数（权值和阈值），训练即告停止。此时经过训练的神经网络即能对类似样本的输入信息，自行处理输出误差最小的经过非线形转换的信息。 2.1 BP模型建立BP网络模型包括其输入输出模型、作用函数模型、误差计算模型和自学习模型。（1）节点输出模型 隐节点输出模型：Oj=f(WijXi-q j) (1) 输出节点输出模型：Yk=f(TjkOj-q k) (2) f-非线形作用函数；q -神经单元阈值。（2）作用函数模型作用函数是反映下层输入对上层节点刺激脉冲强度的函数，又称刺激函数，一般取为(0,1)内连续取值Sigmoid函数： f(x)= 1/(1+e-x) （3） （3）误差计算模型误差计算模型是反映神经网络期望输出与计算输出之间误差大小的函数：Ep=1/2(tpi-Opi)2 tpi- i节点的期望输出值；Opi-i节点计算输出值。 (4)（4）自学习模型神经网络的学习过程，即连接下层节点和上层节点之间的权重拒阵Wij的设定和误差修正过程。BP网络有师学习方式-需要设定期望值和无师学习方式-只需输入模式之分。自学习模型为 Wij(n+1)= hiOj+aWij(n) (5)h-学习因子；i-输出节点i的计算误差；Oj-输出节点j的计算输出；a-动量因子。2.2 BP网络模型的缺陷分析及优化策略（1）学习因子h 的优化采用变步长法根据输出误差大小自动调整学习因子，来减少迭代次数和加快收敛速度。h =h+a(Ep(n)- Ep(n-1)/ Ep(n) a为调整步长，01之间取值 （6）（2）隐层节点数的优化隐层节点数的多少对网络性能的影响较大，当隐节点数太多时，会导致网络学习时间过长，甚至不能收敛；而当隐节点数过小时，网络的容错能力差。利用逐步回归分析 法并进行参数的显著性检验来动态删除一些线形相关的隐节点，节点删除标准：当由该节点出发指向下一层节点的所有权值和阈值均落于死区（通常取0.1、0.05等区间）之中，则该节点可删除。最佳隐节点数L可参考下面公式计算：L=(m+n)1/2+c (7)m-输入节点数；n-输出节点数；c-介于110的常数。（3）输入和输出神经元的确定利用多元回归分析法对神经网络的输入参数进行处理，删除相关性强的输入参数，来减少输入节点数。（4）算法优化由于BP算法采用的是剃度下降法，因而易陷于局部最小并且训练时间较长。用基于生物免疫机制地既能全局搜索又能避免未成熟收敛的免疫遗传算法IGA取代传统BP算法来克服此缺点。2.3 BP神经网络仿真假设备控对象的近似数学模型为： （8）式中，系统a（k）是慢时变的，。神经网络的结构选择2-6-1，学习速度 =0.5和惯性系数=0.05，输入指令信号为：rin(k)=0.5*sin(6*t)。仿真程序：clear all;close all;xite=0.50;alfa=0.05;w2=rand(6,1);w2_1=w2;w2_2=w2;w1=rand(2,6);w1_1=w1;w1_2=w1;dw1=0*w1;x=0,0;u_1=0;y_1=0;I=0,0,0,0,0,0;Iout=0,0,0,0,0,0;FI=0,0,0,0,0,0;ts=0.001;for k=1:1:1000time(k)=k*ts;u(k)=0.50*sin(3*2*pi*k*ts);a(k)=1.2*(1-0.8*exp(-0.1*k);y(k)=a(k)*y_1/(1+y_12)+u_1;for j=1:1:6 I(j)=x*w1(:,j); Iout(j)=1/(1+exp(-I(j);endyn(k)=w2*Iout;e(k)=y(k)-yn(k);w2=w2_1+(xite*e(k)*Iout+alfa*(w2_1-w2_2);for j=1:1:6 FI(j)=exp(-I(j)/(1+exp(-I(j)2;endfor i=1:1:2 for j=1:1:6 dw1(i,j)=e(k)*xite*FI(j)*w2(j)*x(i); endendw1=w1_1+dw1+alfa*(w1_1-w1_2);x(1)=u(k);x(2)=y(k);w1_2=w1_1;w1_1=w1;w2_2=w2_1;w2_1=w2;u_1=u(k);y_1=y(k);endsubplot(2,1,1);plot(time,y,r,time,yn,b);xlabel(times);ylabel(实际输出和仿真输出);grid onsubplot(2,1,2);plot(time,y-yn,r);xlabel(times);ylabel(error);grid on图2-1 神经网络仿真图3动态矩阵控制DMC算法分为三部分:预测模型、滚动优化和反馈校正。3.1 预测模型在DMC算法中，首先需要测定对象单位阶跃响应的采样值,i=1,2,。其中，T为采样周期。对于渐进稳定的对象，阶跃响应在某一，后将趋于平稳，以至(iN)与的误差和量化误差及测量误差有相同的数量级，因而可认为，已近似等于阶跃响应的稳态值。这样，对象的动态信息就可以近似用有限集合加以描述。这个集合的参数构成了DMC的模型参数，向量a=称为模型向量，N称为模型时域长度。虽然阶跃响应是一种非参数模型，但由于线性系统具有比例和叠加性质，故利用这组模型参数，已足以预测在任意输入作用下系统在未来时刻的输出值。在t=kT时刻，假如控制量不再变化时系统在未来N个时刻的输出值为,那么，在控制增量作用后系统的输出可由 （3-1）预测，其中表示在t=kT时刻预测的尚无作用时未来N个时刻的系统输出。表示在t=kT时刻预测的有控制增量作用时未来N个时刻的系统输出。为阶跃响应模型向量，其元素为描述系统动态特性的N个阶跃响应系数。式中，上标表示预测，表示在t=kT时刻预测t=(k+i)T时刻。同样，如果考虑到现在和未来M个时刻控制增量的变化，在t=kT时刻预测在控制增量,作用下系统在未来P个时刻的输出为 （3-2）式中为t=kT时刻预测的无控制增量时未来P个时刻的系统输出。为t=kT时刻预测的有M个控制增量,，时未来P个时刻的系统输出。为从现在起M个时刻的控制增量。称为动态矩阵，其元素为描述系统动态特性的阶跃响应系数。3.2 滚动优化DMC是一种以优化确定控制策略的算法。在采样时刻t=kT的优化性能指标可取为+ （3-3）即通过选择该时刻起M个时刻的控制增量,，使系统在未来P（）个时刻的输出值， 尽可能接近其期望值，。性能指标中的第二项是对控制增量的约束，即不允许控制量的变化过于剧烈。式中，为权系数，P和M分别称为优化时域长度和控制时域长度。显然，在不同时刻，优化性能指标是不同的，但其相对形式却是一致的，都具有类似于(3-1)的形式，所谓“滚动优化”，就是指优化时域随时间不断地向前推移。引入向量和矩阵记号， Q=diag(,，)， R=diag(,，)则优化性能指标式(2-1)可改写为 (3-4)式中，Q, R分别称为误差权矩阵和控制权矩阵。在不考虑输入输出约束的情况下，在t=kT时刻，均为己知，使J(k)取最小的可通过极值必要条件求得 (3-5)这就是t=kT时刻解得的最优控制增量序列。由于这一最优解完全是基于预测模型求得的因而是开环最优解.3.3 反馈校正由于模型误差、弱非线性特性及其它在实际过程中存在的不确定因素，按预测模型式(3-2)得到的开环最优控制规律式(3-5)不一定能导致系统输出紧密地跟随期望值，它也不能顾及对象受到的扰动。为了纠正模型预测与实际的不一致，必须及时地利用过程的误差信息对输出预测值进行修正，而不应等到这M个控制增量都实施后再作校正。为此，在t=kT时刻首先实施中的第一个控制作用 = (3-6) (3-7)其中 (3-8)由于已作用于对象，对系统未来输出的预测便要叠加上产生的影响，即由式(2-1)算出。到下一个采样时刻t=(k+1)T，不是继续实施最优解中的第二个分量，而是检测系统的实际输出y(k + 1)，并与按模型预测算得的该时刻输出，即中的第一个分量进行比较，构成预测误差。 (3-9)这一误差反映了模型中未包含的各种不确定因素，如模型失配、干扰等。由于预测误差的存在，以后各时刻输出值的预侧也应在模型预测的基础上加以校正，这些未来误差的预测，可通过对现时误差e(k+ 1) 加权系数h,(i =1,2,N )得到 （3-10）式中为t=(k+1)T时刻经误差校正后所预测的系统在t=(k+i)T（i=1,，N）时刻的输出.为误差校正向量，其中， =1。 经校正后的的各分量中。除第一项外.其余各项分别是t=(k+1)T时刻在尚无等未来控制增量作用时对输出在t=(k+2)T,(k+N)T时刻的预测值，它们可作为t=(k+1)T时刻的前N一1个分量，即,i=1,，N-1而中的最后一个分量。即t=(k+1)T时刻对i=(k+1+N)T输出的预测，可由来近似，即=，上述关系可用向量形式表示 (3-11)其中为移位矩阵。在t=(k+1)T时刻，有了,就又可以像上面所述t=kT时刻那样进行新的预测优化，整个控制就是在这样推移的过程中滚动进行.由此可以看到，整个动态矩阵控制算法是由预测、控制、校正三部分组成的，该算法结构可用图2-1加以描述。图中粗箭头表示向量流，细箭头表示纯量流。在每一个采样时刻未来P个时刻的期望输出与预测输出所构成的偏差向量按式(3-6)与动态向量点乘。得到该时刻的控制增量。这一控制增量一方面通过数字积分(累加)运算求出控制量u(k)作用于对象;另一方面与阶跃响应向量a相乘，并按式(3-1)计算出在其作用后所预测的系统输出。到了下一个采样时刻，首先测定系统的实际输出y(k+1)，并与原来预测的该时刻的值相比较。按(3-9)算出预测误差e(k+1)。这一误差与校正向量h相乘后，再按式(3-10)校正预测的输出值。由于时间的推移，经校正的预测输出将按式(3-11)移位，并置定为该时刻的预测初值.图中的，表示时移算子，如果把新的时刻重新定义为k时刻，则预测初值的前P个分量将与期望输出一起，参与新时刻控制增量的计算。如此循环，整个过程将反复在线进行。DMC算法是一种增量算法。可以证明，不管有否模型误差，它总能将系统输出调节到期望值而不产生静差。对于作用在对象输入端的阶跃形式的扰动，该算法也总能使系统输出回复到原来的设定状态.3.4 有约束多变量动态矩阵控制及其线性化不失一般性, 假设阶跃响应的截断时域长度均为N , 控制时域长度均为M , 预测时域长度均为P。基于阶跃响应的预测模型为 (3-12)其中为k时刻第j个输入对第i个输出的贡献，为第i个输出对第j个输入的第s个阶跃响应系数。则第i个模型输出为 （3-13）令，考虑到反馈校正作用，把控制量分解为当前及将来未知部分和过去已知部分。可得k时刻输出的P步预测为： （3-14）其中，和分别为和阶矩阵，其元素由决定，为维校正列向量，而将化成全量形式，有 （3-15），将式（3-15）代入（3-14）得 （3-16）令，则式（2-16）成为 （3-17）设第i个过程的P步参数输出为为跟踪参考轨迹，应使为使上式性能指标线性化，将其改写成 （3-18）其中，和表示相应维数的正、负偏差变量列向量，对此二向量求其最小指对其范数求最小。令结合式(2-18)和（2-17）可得： （3-19）对控制施加幅值约束，不失一般性，假设控制幅值的下限为零，即 （3-20）其中，为第j个控制的最大幅值，。如果控制幅值的下限约束为负，则可通过坐标或利用线性规划变量非负化的方法使其为零，以满足单纯形法迭代寻优的要求。对式（3-20）引入相应维数的正、负偏差变量列向量，使不希望偏差最小化，则式（3-20）为： （3-21）对控制变化速率施加上、下限约束, 即 （3-22）其中，为第j 个控制变化速率的上、下限。同理, 引入相应维数的正、负偏差变量列向量, 最小化不希望偏差, 则式(2-22) 可表示为 （3-23）3.5 动态矩阵控制仿真设备控对象的近似数学模型为： 式中，系统a（k）是慢时变的，输入指令信号为：rin(k)=10。 仿真程序：clear all;close all;u_1=0;y_1=0;for k=1:1:1000u(k)=1;a(k)=1.2*(1-0.8*exp(-0.1*k);yg(k)=a(k)*y_1/(1+y_12)+u_1;y0(k)=yg(k);y_1=yg(k);u_1=u(k);endTs=5;N=200;P=25;m=20;Sv=10;yg_1=0;A=zeros(P,m);a=zeros(N,1);for i=1:N a(i)=y0(i);endfor i=1:P for j=1:m if i-j+10 A(i,j)=a(i-j+1); end endendK=inv(A*A+eye(m)*A;y=zeros(N,1);u=zeros(N,1);e=zeros(N,1);A0=zeros(P,N-1);for i=1:P for j=N-2:-1:1 if N-j+1+i-1=N A0(i,j)=a(N-j+1+i-1)-a(N-j+i-1); else A0(i,j)=0; end end A0(i,N-1)=a(i+1);endfor k=2:N Uk_1=zeros(N-1,1); for i=1:N-1 if k-N+i=0 Uk_1(i)=0; else Uk_1(i)=u(k-N+i); end end Y0=A0*Uk_1;a(k)=1.2*(1-0.8*exp(-0.1*k);yg(k)=a(k)*y_1/(1+y_12)+u_1;u_1=u(k); yg_1=yg(k); e(k)=y(k)-yg(k); Ysk=zeros(P,1); for i=1:P Ysk(i)=Sv; end Ek=zeros(P,1); for i=1:P Ek(i)=0.1*i*e(k); end dertu=K*(Ysk-Y0-Ek); for i=1:m if k+i-1=N u(k+i-1)=u(k+i-1-1)+dertu(i); end end temp=0; for j=1:N-1 if k-j=0 temp; else if k-j-1=0 temp=temp+a(j)*u(k-j); else temp=temp+a(j)*(u(k-j)-u(k-j-1); end end end if k-N0 A(i,j)=a(i-j+1); end endendK=inv(A*A+eye(m)*A;y=zeros(N,1);u=zeros(N,1);e=zeros(N,1);A0=zeros(P,N-1);for i=1:P for j=N-2:-1:1 if N-j+1+i-1=N A0(i,j)=a(N-j+1+i-1)-a(N-j+i-1); else A0(i,j)=0; end end A0(i,N-1)=a(i+1);endfor k=2:N Uk_1=zeros(N-1,1); for i=1:N-1 if k-N+i=0 Uk_1(i)=0; else Uk_1(i)=u(k-N+i); end end Y0=A0*Uk_1;a(k)=1.2*(1-0.8*exp(-0.1*k);yg(k)=a(k)*y_1/(1+y_12)+u_1;u_1=u(k); yg_1=yg(k); e(k)=y(k)-yg(k); Ysk=zeros(P,1); for i=1:P Ysk(i)=Sv; end Ek=zeros(P,1); for i=1:P Ek(i)=0.1*i*e(k); end dertu=K*(Ysk-Y0-Ek); for i=1:m if k+i-1=N u(k+i-1)=u(k+i-1-1)+dertu(i); end end temp=0; for j=1:N-1 if k-j=0 temp; else if k-j-1=0 temp=temp+a(j)*u(k-j); else temp=temp+a(j)*(u(k-j)-u(k-j-1); end end end if k-N=0 y(k)=temp+e(N); else y(k)=temp+a(N)*u(k-N)+e(N); endendt=Ts.*(1:N);figure(3);subplot(2,1,1);plot(t,y,b,t,yg,r);legend(预测输出曲线);title(输出曲线);xlabel(Time)ylabel(振幅)grid onsubplot(2,1,2);plot(t,u ,g);legend(控制作用u)title(控制作用);xlabel(Time)ylabel(振幅)grid onxite=0.50;alfa=0.05;w2=rand(10,1);w2_1=w2;w2_2=w2;w1=rand(6,10);w1_1=w1;w1_2=w1;dw1=0*w1;u_1=0;y_1=0;I=zeros(10,1);Iout=zeros(10,1);FI=zeros(10,1);x=0,0,0,0,0,0;for k=3:N x(1)=yg(k); x(2)=yg(k-1); x(3)=yg(k-2); x(4)=u(k); x(5)=u(k-1); x(6)=u(k-2);for j=1:1:10 I(j)=x*w1(:,j); Iout(j)=1/(1+exp(-I(j);endes(k)=w2*Iout;e1(k)=es(k)-e(k);w2=w2_1-(xite*0.1*e1(k)*Iout+alfa*(w2_1-w2_2);for j=1:1:10 FI(j)=exp(-I(j)/(1+exp(-I(j)2;endfor i=1:1:6 for j=1:1:10 dw1(i,j)=e1(k)*xite*FI(j)*w2(j)*x(i); endendw1=w1_1+dw1+alfa*(w1_1-w1_2);w1_2=w1_1;w1_1=w1;w2_2=w2_1;w2_1=w2;endt=Ts.*(1:N);figure(2);plot(t,es,r,t,e,b);xlabel(times);ylabel(es and e);grid onfor k=3:N Uk_1=zeros(N-1,1); for i=1:N-1 if k-N+i=0 Uk_1(i)=0; else Uk_1(i)=u(k-N+i); end end Y0=A0*Uk_