多入多出非线性预测控制

多入多出非线性预测控制，可拿来做研究和学习。

clear all
load wentaishuju.mat;
load ARX.mat;
load BPmoxing.mat;
%-----------参数设定------------------%
nu=2;    %系统的输入维数
ny=2; 
P=5;
M=3;
c=0.5;
d=0;
P>d;
（P-d）>=M;  %所给出的初始值一定要满足的条件
miu=0.4;   %  miu是调节模型参数的遗忘因子在0～1范围内调节
a=0.9;
I=eye（ny）;
dus=[0.04;0.04];   %用来求解增益Kc和Kn的输入量的小的变化量
x0=[1;1;1;1;1;1];
s1=[40 40 40 40 40 40 40 40 40 40];
s2=[10 10 10 10 10 10];
Q=diag（s1）;
R=diag（s2）;
%H=[0.1*I;0.1*I;0.1*I;0.1*I;0.1*I];% 暂定为0.1*I
%H1=ones（51）;
%H2=zeros（51）;
H=[I;I;I;I;I];
%H=[H1 H2;H2 H1];  %修改了H后 与原文不同
%-----------参考轨迹------------------%
for i=1:3;
    y1sp（i）=0;
end
for i=4:100;
    y1sp（i）=2;
end
for i=101:200;
    y1sp（i）=2;
end
for i=201:250;
    y1sp（i）=3; 
end
for i=251:300;
    y1sp（i）=3;
end
for i=301:400;
    y1sp（i）=2;
end


for i=1:20;
    y2sp（i）=0;
end
for i=21:150;
    y2sp（i）=3;
end
for i=151:250;
    y2sp（i）=3;
end
for i=251:320;
    y2sp（i）=2; 
end
for i=321:400;
    y2sp（i）=3;
end
ysp=[y1sp;y2sp];  %两输出设定值组合矩阵
for i=1:ny;
%Vr（ik）=vr（k+i）;
C1（i1）=c;
%C1（i1）=1-c^i;
end
C=diag（C1）;



%-----------初始模型的参数Kc----------------%
%{
for i=1:NO;                                   %
      if （yr（1）>=ys（i））&&（yr（1）<=ys（i+1））;
          ysg=ys（i+1）;
          ysl=ys（i）;
          break
      end
end
 Kc1=（ysg-ysl）/dumin;
%theta（:1）=[a11;a12;b1*Kc1;b2*Kc1]; %Kc1是us（1）对应的动态增益
for i=1:NO;
      if （yr（2）>=ys（i））&&（yr（2）<=ys（i+1））;
         ysg=ys（i+1）;
         ysl=ys（i）;
          break
      end
end
 Kc2=（ysg-ysl）/dumin;
 %theta（:2）=[a11;a12;b1*Kc2;b2*Kc2]; %Kc2是us（2）对应的动态增益
 
 I=[1 0 0 0;0 1 0 0;0 0 1 0;0 0 0 1];
P1=（a^2）*I;

%}
%---------------预测循环开始-------------%
for i=1:2;
yc（1:2i）=0.0001;
%v（i）=0.0001;
%u（i）=0.0001;
%x（i）=0.0001;
end
y（:1）=[0;0];
u（:1）=[1;0];
y（:2）=[0;0];
y（:3）=[0;0];
u（:2）=[1;0];
du（:1）=[0;0];
du（:2）=[0;0];

step=300;
for k=2:step;
    %{
        u11=u（1k-1）;
        u12=u（2k-1）;
        [tx1]=ode45（‘myssFun‘[0 500][40 0.17]）;%初始化状态变量
        [hl]=size（x1）;
        x1（h1）=abs（x1（h1））;
        x1（h2）=abs（x1（h2））;
        y1（k）=x1（h1）;       % 单输出数据
        y2（k）=x1（h2）;
        y（:k）=[y1（k）;y2（k）]; % 混合输出数据的对应矩阵
   
    uz=u（1k-1）;
    u（1k-1）=u（2k-1）;
    u（2k-1）=uz;
     %}
    y（:k）=sim（network1u（:k-1））;
dy（:k）=y（:k）-y（:k-1）;
%Yr（:k）=C*y（:k）+C1*ysp（:k+d）;
yr（:k）=y（:k）;  %第k步的参考轨迹是第k步的实际输出值
for i=1:P;
    yr（:k+i）=C*yr（:k+i-1）+（I-C）*ysp（:k）;
    [hl]=size（yr（:k+i））;
    Yr（（i*h-h+1）:i*hk）=yr（:k+i）;  %Yr是一个（ny*P）*1的列向量
end

%for i=2:NO-1;
    %  if （y（k）>=ys（i））&&（y（k）<=ys（i+1））;
    %     ysg=ys（i+1）;
    %     ysl=ys（i）;
        
    %      break
     % end
%end
%usc=YU（i2）;
% Kc=（ysg-ysl）/dumin;   %求得滚动变化的Kc
% for i=2:NO-1;
  %    if （yr（k+1）>=ys（i））&&（yr（k+1）<=ys（i+1））;
  %       ysg=ys（i+1）;
 %        ysl=ys（i）;
 %         break
  %    end
% end
%usn=YU（i2）;
% Kn=（ysg-ysl）/dumin;   %求得滚动变化的Kn


%----------------

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     文件      148266  2010-08-26 14:59  BPmoxing.mat
     文件       84013  2013-12-24 15:57  wentaishuju.mat
     文件       15520  2010-12-01 11:38  zuhemoxingcontrol.m
     文件         448  2010-11-06 14:59  ARX.mat