dmc.m

来自「很好的预测控制DMC程序」· M 代码 · 共 110 行

clear all;

%设置初始值
N=20;           
q=0.5;          
r=0.4;          
M=2;            
P=10;           
ysp=1;          
rank=5;

%计算动态矩阵A
for k=1:N
    u(k)=1;
end

a=[];
for k=1:rank;           
    a=[a,0];
end

for k=1:N               
    atemp=0.2713*u(1)+0.8351*a(k+4);
    a=[a,atemp];
end

Atemp=[];
A=[];
for k=1:M
    Atemp=[Atemp;zeros(1:k-1),a(1:N+1-k)];
    A=Atemp';
end

%计算控制器参数d
Q=eye(P)*q;
R=eye(M)*r;
d=[];
A1=A(1:P,1:M);
d=(inv(A1'*Q*A1+R)*A1'*Q)';

%设置反馈系数h，参考轨迹系数alpa
h=1.5;
alfa=0.8;

%初始化DMC
Alfa=[];
for j=1:P
    Alfa=[Alfa,alfa^j];
end

H=[];                           
for k=1:N
    H=[H;h];
end

Y=[0];
Y0=[];
Y1=[];
for k=1:P                  
    Y0=[Y0;0];   
end

Yn0=[];
for k=1:N
    Yn0=[Yn0;0];
end

Yn1=[];                    
U=[];
Delta_U=[];

%模拟控制输入和输出,检验控制效果
Utemp=[0;0];
U1=[];

for k=1:100
    W=(1-Alfa')*ysp+Alfa'*Y(k); 
    
    delta_U=d'*[W-Y0];          
    Delta_U=[Delta_U,delta_U];  
    Utemp=Utemp+delta_U;        
    U=[U,Utemp];

    if k<=rank
        y=0;
    else
        y=10.2713*U(1,k)+0.8351*Y(k);
    end
    Y=[Y,y];
    
    Y1=Y0+A(1:P,1:M)*Delta_U(:,k);
    error=Y(k+1)-Y1(1,1);
    
    Yn1=Yn0+A(1:N,1:M)*Delta_U(:,k);
    
    Ycor=[Yn1+H*error];

    Yn0=[Ycor(2:N);Ycor(N)];
    Y0=[Yn0(1:P)];
end
Y=Y(2:k+1);
t=1:100;
figure(1)
plot(t,Y)
xlabel('N');
ylabel('Y');
figure(2)
plot(t,U)
xlabel('N');
ylabel('U');