nnpc神经网络预测控制matlab代码（单步预测）

clear all; close all; a=[1 -0.7205]; B=[0.6849]; %a=[1.0000 -0.9677]; %B=[0.0790]; %c=1; %d=2; %对象参数 %a=[1 -1.4785 0.5326]; %B=[0.018 0.0146]; %c=1; global d he; he=0; d=5 %对象参数 na=length(a)-1; B=[zeros(1,d-1) B]; nb=length(B)-1; %na、nb为多项式A、B阶次（因d!=1，对b添0） aa=conv(a,[1 -1]); naa=na+1; %aa的阶次 global N N1=d; N=2; Nu=5; %最小输出长度、预测长度、控制长度 gamma=1*eye(Nu); alpha=0.6; %控制加权矩阵、输出柔化系数 L=200; %控制步数 r=2*rand(L+d+10)-1; uk=zeros(nb+1,1); %输入初值：uk(i)表示u(k-i) duk=zeros(d+nb,1); %控制增量初值 yk=zeros(na,1); %输出初值 global W W=5*[ones(L/4,1);-ones(L/4,1);ones(L/4,1);-ones(L/4+d+10,1)]; %设定值 %w=500*[zeros(L/2,1);ones(L/2+d,1)]; %W=5*[ones(L+d+10,1)]; % 100 90 80 70 60 50 40 30 20 10 0 -10 -20 -30 -40 -50 -60 -70 -80 -90 -100; global c c=[100 90 80 70 60 50 40 30 20 10 0 -10 -20 -30 -40 -50 -60 -70 -80 -90 -100; 100 90 80 70 60 50 40 30 20 10 0 -10 -20 -30 -40 -50 -60 -70 -80 -90 -100; ]; global b b=30; global YY U1 YY=zeros(d+N-1,1); U1=zeros(d+N-1,1); %100 90 80 70 60 50 40 30 20 10 0 -10 -20 -30 -40 -50 -60 -70 -80 -90 -100; %w=rands(1,21); global w h2; h3=zeros(21,1); w=rands(1,21); w=[-0.348208037112170,0.725896829339912,0.674097916559497,0.775642844135692,2.14502719673641,4.74999936740020,6.63942397279853,9.47410336309908,9.75015806080596,6.48956982303472,0.323220101083847,-8.93142400002481,-11.3064379976673,-11.5672257732651,-7.76399162247795,-3.95555069150801,-1.75448273611684,0.0413093244632264,0.608958863206572,0.327970649388938,-0.506881799500054] %w=[0.692993678522641,2.94563273762520,4.43498568339060,8.09305084548958,11.6517954193711,13.0333789659908,13.8353390446874,13.0911663950922,9.25901051541135,5.58235946384544,-0.186782503486349,-2.77268159992267,-6.24000640700897,-8.43690742716629,-11.0962617926567,-12.6955888270415,-12.0449163426565,-9.98110615284044,-8.72965557971494,-5.25554145784800,-3.70208899789528] %w=[0.519625597882877,0.248628895865967,1.40206845759801,3.55795723573982,5.04493899178630,6.19033613737719,8.25373934131726,8.83586685420943,8.15411652270847,4.21640651086086,1.09133106896299,-2.94967096645171,-8.67061976104969,-9.80763252883009,-9.64454886334168,-8.94845653005136,-6.39097494188769,-4.82428374789007,-2.59748679596066,-2.15903688736879,-0.780646358957199] xite=0.1; alfa=0.05; h=[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]'; %b=(7/10)*[50,50,50,50,50,50,50,50,50,50,50,50,50,50,50,50,50,50,50,50,50]; c1=c;c2=c1; b1=b; b2=b1; w1=w; w2=w1; Y=0; yr=0; e1=0; y=zeros(L,1); h2=zeros(21,1); global ukd ykd k for k=1:1:L time(k)=k; y(k)=-a(2:naa)*yk+B*uk(1:nb+1); %采集输出数据 % y(k)=sin(0.1*k); %模型辨识 x=[yk',uk(d:nb+1)']; yout=y(k); for j=1:1:21 h(j)=exp(-norm(x'-c(:,j))^2/(2*b^2)); end ymout=w*h; ec(k)=ymout; dw=0*w; for j=1:1:21 for i=1:1:2 c(i,j)=c(i,j)+xite*(yout-ymout)*w(j)*exp(-norm(x'-c(:,j))^2/(2*b^2))*((x(i)-c(i,j))/(b^2))+alfa*(c1(j)-c2(j)); end end b1=b; b2=b1; c1=c;c2=c1; for j=1:1:21 %Only weight value update dw(j)=xite*(yout-ymout)*h(j); end w=w1+dw+alfa*(w1-w2); w2=w1;w1=w; %滚动优化： du=uk(1); dtu=0; yc=y(k)-ymout; yout=ymout+0.1*yc; ykd=yk; ukd=uk; % du=1; for i=1:1:d+N xd=[ykd',ukd(d:nb+1)']; for j=1:1:21 h2(j)=exp(-norm(xd'-c(:,j))^2/(2*b^2)); end yd=w*h2; YY(i)=yd; U1(i)=du; for i=nb+1:-1:2 ukd(i)=ukd(i-1) end ukd(1)=du; % yc=y(k)-ymout; % yout=ymout+0.1*yc; for i=na:-1:2 ykd(i)=ykd(i-1) end ykd(1)=yd; end %xd=[ykd',ukd(d:nb+1)']; % ukd(5)=du; ukd(d)=du; %xd=[ykd',ukd(d)']; %for j=1:1:21 % h3(j)=exp(-norm(xd'-c(:,j))^2/(2*b^2)); %end yd=w*h2; f = @(du)fun(du,YY); rng default nvars = 1; syms du; x0 = [1]; % du=2 [du,fval]=fminunc(f,x0); % [du,fval]=particleswarm(f,2); % du=20*randn(1,1) %fval %he=0; %参考轨迹 %{ for i=1:1:d xd=[ykd',ukd(d:nb+1)']; for j=1:1:21 h2(j)=exp(-norm(xd'-c(:,j))^2/(2*b^2)); end yd=w*h2'; for i=nb+1:-1:2 ukd(i)=ukd(i-1) end ukd(1)=du; for i=na:-1:2 ykd(i)=ykd(i-1) end ykd(1)=yd; end %} %{ yr=0*yr+(1-0)*W(k+d);%期望输出柔画 kp=1; ki=0.05; kd=0.05; yr=0*yr+(1-0)*W(k+d);%期望输出柔画 e=W(k+d)-ykd(1); % e=W(k)-y(k); ed=e-e1; ei=e+e1; e1=e; du=e*kp+ei*ki+ed*kd; U(k)=du; %} %{ for m=1:1:200 x=[ykd',ukd(d:nb+1)']; dtu=0; for i=1:1:21 h1(i)=exp(-norm(x'-c(:,i))^2/(2*b^2)); dtu=dtu+w(i)*h1(i)*((c(2,i)-x(2)))/(b^2); end du=du+0.05*(yr-w*h1')*dtu; ukd(d)=du; % ykd(1)=w*h1'; end %} % if(du<-20) % du=-20; % end % if(du>20) % du=20; % end for i=nb+1:-1:2 uk(i)=uk(i-1) end uk(1)=du; % uk(1)=20*r(k); % uu(k)=uk(1)+5; for i=na:-1:2 yk(i)=yk(i-1); end % yk(2)=yk(1); yk(1)=y(k); end plot(time,y,'r',time,W(1:L),time,ec,'y'); xlabel('k'); ylabel('y(k)');