串级DMC+PID控制

clear all;
close all;
clc;

ts=1;
sys1=tf（1[4.7871 1]）;        %比例系数？？
sys1=sys1*exp（-10）;
dsys1=c2d（sys1ts‘z‘）;
[num1den1]=tfdata（dsys1‘v‘）; %上水箱

u_2=0;
u1_1=0;
y2_1=0;
y1_1=0.1;
ei=0;
x=[000]‘;
error_1=0;
kp=80;ki=2;kd=0; 
% p=13;       %  p=17m=3.%外环
% m=2;
ywt=[];
uwt=1;
N=500;

% yr=50;

for s=1:1:100 

    time（s）=s*ts; 
    %Linear model
    G=tf（[0.1426][536 46.3024 1]）;
    G=G*y1_1;
    P=50;M=20;
    A=zeros（PM）;%动态矩阵
    a=zeros（N1）;

    [y0t0]=step（G0:1:500）;

    for i=1:N
        a（i）=y0（i）;
    end
    for i=1:P
        for j=1:M
            if i-j+1>0
                A（ij）=a（i-j+1）; %构造矩阵A
            end
        end
    end
    %初始化向量ys，yue和矩阵A0
    ys=ones（N1）;
    y=zeros（N1）;
    u=zeros（N1）;
    e=zeros（N1）;
    A0=zeros（PN-1）;
    for i=1:P
        for j=N-2:-1:1
            if N-j+1+i-1<=N
                A0（ij）=a（N-j+1+i-1）-a（N-j+i-1）;%构造矩阵A0
            else
                A0（ij）=0;
            end
        end
        A0（iN-1）=a（i+1）;
    end
    %DMC程序
    for k=2:N
        if k>=N/2&&k<=N/2+3
            yr=60;
        else
            yr=50;
        end
        Uk_1=zeros（N-11）;
        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;
        e（k）=y（k-1）-Y0（1）;
        Yr=zeros（P1）;
        for i=1:P
            Yr（i）=yr;
        end
        Ek=zeros（P1）;
        for i=1:P
            Ek（i）=e（k）;
        end
        delta_u=inv（A‘*A+eye（M））*A‘*（Yr-Y0-Ek）; %控制增量的计算
        for i=1:M
            if k+i-1<=N
                u（k+i-1）=u（k+i-1-1）+delta_u（i）;  %控制律的计算
            end
        end
        temp=0;%设置在k-j-1时刻以前的控制律
        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）;
        end
    end
 
    error（s）=u（s）-y1_1;
    x（1）=error（s）;                 %Calculating P
    x（2）=（error（s）-error_1）/ts;  %Calculating D
    x（3）=x（3）+error（s）*ts;        %Calculating I

    u1（s）=kp*x（1）+kd*x（2）+ki*x（3）; 
    y1（s）=-den1（2）*y1_1+num1（2）*u1（s）;  %内环
    % 
    error_1=error（s）;
    % u1_1=u1（k）;
    y1_1=y1（s）;
    % u_2=u（k）;
    y2_1=y（s）;

end
 t=（1:N）;
% subplot（211）;
plot（ty‘r‘）
hold on;
% title（‘PID串级PID和预测PID输出‘）;
xlabel（‘t‘）
ylabel（‘y‘）
%grid on
% subplot（212）;
% plot（tu‘r‘）;
% title（‘控制作用‘）;
% xlabel（‘t‘）
% ylabel（‘u‘）
%grid on
ts=1;
sys1=tf（0.1426[4.78711]）;

dsys1=c2d（sys1ts‘z‘）;
[num1den1]=tfdata（dsys1‘v‘）;

sys2=tf（1[23.04941]）;
% sys2=sys2*exp（-10）;
dsys2=c2d（sys2ts‘z‘）;
[num2den2]=tfdata（dsy