MPC_DMC.rar_MPC matlab_MPC预测控制_matlab mpc_预测控制

clc; clear all; STEPS=240; %%预测步长 P=30; %%控制计算步长 M=10; g=zeros(1,STEPS); u=ones(1,STEPS); for i=1:STEPS if i>1 g(i)=0.9943*g(i-1)+2.747*u(i-1); end end % g(1:30) =[ 0; 0;0.2713;0.4979;0.6871;0.8451;0.9770;1.0872;1.1792;1.2561;1.3202;1.3738; % 1.4186;1.4560;1.4872;1.5132;1.5350;1.5532;1.5684;1.5810;1.5916;1.6005;1.6079; % 1.6140;1.6192;1.6235;1.6271;1.6301;1.6326;1.6346]; % g(1:40)=[9.4505e-011 1.2348e-006 0.00018178 0.0038545 0.02624 0.083386 0.15891 0.23011 0.30511 0.38655 0.46194 0.53078 0.59576 0.654 0.70507 0.75015 0.78923 0.82269 0.85123 0.87537 0.89566 0.91261 0.92671 0.93837 0.94799 0.95589 0.96235 0.96761 0.9719 0.97537 0.97818 0.98044 0.98227 0.98373 0.98491 0.98585 0.9866 0.9872 0.98767 0.98805]; % for i=41:1:STEPS % g(i)=g(40); % end G=zeros(P,M); for i=1:M G(i:P,i)=(g(1:P-i+1))'; end %%控制权系数 Q=10; R=2; d=inv(G'*Q*eye(P)*G+R*eye(M))*G'; %w0=ones(P,1); w=[]; for i=50:1:150 w(i)=1; end for i=151:STEPS+P w(i)=2; end wtemp_2=zeros(P,1); yout=zeros(1,STEPS); ytemp=0; ucon=[]; de_u=[]; for i=1:1:STEPS ucon(i)=0; de_u(i)=0; end upred=zeros(M,1); utemp=0; %%当前使用的控制量 ystates=[]; ypred2=[]; dr_flag=[]; delay=[]; delay1=[]; dr_order=[]; ylast=0; for i=1:1:STEPS+5 ystates(i)=0; dr_flag(i)=0; ypred2(i)=0; delay(i)=0; dr_order(i)=0; end max_order=dr_order(1); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for counter=1:1:STEPS delay(counter)=rand(1,1)*4+1; delay1(counter)=fix(delay(counter)); for i=1:1:counter-1 if dr_order(i)>max_order max_order=dr_order(i); end end % 利用t-1时刻及以前时刻控制器接收到的最新测量值的时间标志,与t时刻接收到的数据时间标志进行比较. % 若有更新的数据到达,则控制器接收端数据更新为此数据 % 若没有新的数据到达，则控制器接收端数据更新为模型预测的当前值 if dr_order(counter)>=max_order && dr_flag(counter)==1 %% dr_flag(counter)==1表示当前时刻有数据到达控制器侧 max_order=max(dr_order(counter)); ylast=ystates(counter); else ylast=ypred2; %% 对应数据包丢失的情况以及数据包乱序的情况 end for i=1:1:P wtemp_2(i)=w(counter+i); end %下面计算自由响应向量 fprev=zeros(P,1); for i=1:1:P for j=1:1:200 if (counter-j)>0 temp=(g(i+delay1(counter)+j)-g(j))*de_u(counter-j); fprev(i)=fprev(i)+temp; end end end if (counter-1)>0 fprev=fprev+(ylast*ones(P,1)); end % 下面计算当前时刻的最优控制量 upred=d*(wtemp_2-fprev); % 仅使用第一个控制量 de_u(counter)=upred(1); if (counter-1>0) utemp=upred(1)+ucon(counter-1); else utemp=upred(1); end ucon(counter)=utemp; %计算在最优控制量作用下的实际输出 if (counter-1)>0 ytemp=0.9943*ylast+2.747*ucon(counter-1); yout(counter)=ytemp; end ystates(counter+delay1(counter))=ytemp; dr_flag(counter+delay1(counter))=1; dr_order(counter+delay1(counter))=counter; % 用计算出的最优控制量upred计算预测输出 ypred=G*upred+fprev; ypred2=ypred(1); end subplot(3,1,1) stairs(ucon) ylabel('u') subplot(3,1,2) plot(yout) ylabel('y') hold on plot(w,':') subplot(3,1,3) plot(delay) ylabel('delay')