模型预测控制（MPC）MATLAB实现-01

%%%%%% 模型预测控制MPC的MATLAB实现

clear;clc;

%%%% 系统模型辨识函数
%ex8_1 对单输入单输出系统的传递函数进行脉冲响应模型辨识
num=[1 1];
den=[1 3 6];
h=tf(num,den); %转化为传递函数模型

[u,t]=gensig('pulse',2,10,0.1); %获得指定脉冲信号u/x
x=u;
%求解LTI对象的单位脉冲响应y
t=0:0.1:10;
[y,xl,tl]=lsim(h,x,t); %获得任意输入函数的输出响应

%对输入脉冲信号x进行归一化处理
[ax,mx,stdx]=autosc(x); 
mx=0;
sx=scal(x,mx,stdx);

n=35; %脉冲响应模型的个数
[xreg,yreg]=wrtreg(sx,y,n); %生成用于线性回归计算的输入输出数据矩阵

%基于多变量最小二乘法的脉冲响应模型辨识函数mlr()。
ninput=1;
plotopt=2;
[theta,yres]=mlr(xreg,yreg,ninput,plotopt);

%ex8_2 对双输入单输出系统的传递函数进行脉冲响应模型辨识
num1=5.72;den1=[60 1];g1=poly2tfd(num1,den1,0,14);
num2=1.52;den2=[25 1];g2=poly2tfd(num2,den2,0,15); % 传递函数转换为MPC传递函数
mod=tfd2mod(7,1,g1,g2); % MPC传递函数转换为MPC状态空间模型
[A,B,C,D]=mod2ss(mod); % MPC状态空间模型转换为状态空间模型

%将通用状态空间模型转换为LTI对象
sys=ss(A,B,C,D);h=tf(sys);

[u,t]=gensig('pulse',8,10,0.1); 
x=[u u];
%求解LTI对象的单位脉冲响应y
t=0:0.1:10;
[y,xl,tl]=lsim(h,x,t); %获得任意输入函数的输出响应

%对输入脉冲信号x进行归一化处理
[ax,mx,stdx]=autosc(x); 
mx=[0 0];
sx=scal(x,mx,stdx);

n=35; %脉冲响应模型的个数
[xreg,yreg]=wrtreg(sx,y,n); %生成用于线性回归计算的输入输出数据矩阵

%基于多变量最小二乘法的脉冲响应模型辨识函数mlr()。
ninput=2;
plotopt=2;
[theta,yres]=mlr(xreg,yreg,ninput,plotopt);


%ex8_3 单输入单输出传递函数利用部分最小二乘法进行脉冲响应模型辨识
num=[1 1];
den=[1 3 6];
h=tf(num,den); %转化为MPC传递函数模型

[u,t]=gensig('pulse',2,10,0.1); %获得指定脉冲信号u/x
x=[u];

t=0:0.1:10;
[y,xl,tl]=lsim(h,x,t); %获得任意输入函数的输出响应

n=30; %脉冲响应模型的个数
[xreg,yreg]=wrtreg(sx,y,n); %生成用于线性回归计算的输入输出数据矩阵

%基于部分最小二乘法（PLS）的脉冲响应模型辨识
ninput=2;
lv=10;
plotopt=2;
theta=plsr(xreg,yreg,ninput,lv,plotopt);

%ex8_4 脉冲响应模型转换为阶跃响应模型 
num=[1 1];den=[1 3 6];h=tf(num,den); %传递函数模型转换为MPC传递函数模型
[u,t]=gensig('pulse',2,10,0.1); x=[u]; %获得脉冲信号x
%求解LTI对象的单位脉冲响应y
t=0:0.1:10;
[y,xl,tl]=lsim(h,x,t); %获得任意输入函数的输出响应
%对输入脉冲信号x进行归一化处理
[ax,mx,stdx]=autosc(x); mx=[0];sx=scal(x,mx,stdx);

n=35; %脉冲响应模型的个数
[xreg,yreg]=wrtreg(sx,y,n); %生成用于线性回归计算的输入输出数据矩阵

%基于多变量最小二乘法的脉冲响应模型辨识
ninput=1;
plotopt=2;
[theta,yres]=mlr(xreg,yreg,ninput,plotopt);
%脉冲响应模型转换为阶跃响应模型
theta=scal(theta,mx,stdx);nout=1;delt=1;
model=imp2step(delt,nout,theta);
plotstep(model); %系统阶跃响应函数


%%%% 系统模型建立与转换函数
%ex8_5 传递函数转换为MPC状态空间模型
num=[1 3 1];den=[1 2 2 1];
[A,B,C,D]=tf2ss(num,den); % 传递函数转换为状态空间方程
pmod=ss2mod(A,B,C,D); % 状态空间方程转换为MPC状态空间方程

[A,B,C,D]=mod2ss(pmod); % MPC状态空间方程转换为状态空间方程

%ex8_6 传递函数模型转换为MPC传递函数模型
num=[1 1];den=[1 4 4];
g=poly2tfd(num,den,0,0.5); % 传递函数转换为MPC传递函数

%ex8_7 MPC传递函数转换为MPC状态空间模型
num=[1 1];den=[1 3 6];
g=poly2tfd(num,den,0,0); % 传递函数转换为MPC传递函数
modl=tfd2mod(0.1,1,g); % MPC传递函数转换为MPC状态空间模型

%ex8_8 传递函数模型转换为阶跃响应模型
num=[1 2];den=[1 3 1];
tfl=poly2tfd(num,den,0,0); % 传递函数转换为MPC传递函数
plant=tfd2step(5,0.1,1,tfl); % MPC传递函数转换为阶跃响应模型
plotstep(plant); %系统阶跃响应曲线


%ex8_9 多变量系统传递函数模型转换为阶跃响应模型
tf1=poly2tfd(2,[1,1],0,0);tf2=poly2tfd(3,[1,2],0,0);
tf3=poly2tfd(1,[1,4],0,0);tf4=poly2tfd(1,[1,1],0,0); % 传递函数模型转换为MPC传递函数模型
plant=tfd2step(5,0.1,2,tf1,tf2,tf3,tf4); % MPC传递函数模型转为阶跃响应模型
plotstep(plant); % 逐步执行

%ex8_10 状态空间模型转换为阶跃响应模型
A=[-2 -1;1 0];B=[1;0];C=[1 2];D=0;
plant=ss2step(A,B,C,D,5,0,0.2,1); % 状态空间模型转换为阶跃响应模型,截断时间5s，阶跃响应的采样周期为0.2
plotstep(plant); % 系统的阶跃响应曲线

%ex8_11 改变MPC状态空间模型的采样周期
num=[1 1];den=[1 3 6];
g=poly2tfd(num,den,0,0); % 传递函数转换为MPC传递函数

mod1=tfd2mod(0.1,1,g); % 采样周期，输出个数，MPC传递函数

% plant1=mod2step(mod1,5,0.1);
% plotstep(plant1)

mod2=mod2mod(mod1,0.5); % 离散系统原有的MPC状态空间方程 指定系统新的采样周期
plant2=mod2step(mod2,5,0.1);
plotstep(plant2)


%%%% 基于阶跃响应模型的控制器设计与仿真函数
%ex8_12 已知系统传递函数，对其进行输入/输出有约束的模型预测控制器设计仿真
%将传递函数模型转换为MPC传递函数模型
g11=poly2tfd(12.8,[16.7 1],0,1);   g21=poly2tfd(6.6,[10.9 1],0,7);
g12=poly2tfd(-18.9,[21.0 1],0,3);  g22=poly2tfd(-19.4,[14.4 1],0,3);

% g=polytfd(num,den,delt,delay)
% delt 为采样周期，，对于连续系统，该参数为0,；delay 为系统纯延时，对于离散系统，纯延时是采样周期的整数倍；

tfinal=90;
delt=3;
ny=2;
model=tfd2step(tfinal,delt,ny,g11,g21,g12,g22);

%% [y,u,ym]=cmpc(plant,model,ywt,uwt,M,P,tend,r,ulim,ylim，tfilter,dplant,dmodel,dstep);
%  plant为开环对象的实际阶跃响应模型；model为辨识得到的开环对象阶跃响应模型；
%  ywt为二次型性能指标输出误差加权矩阵Q；uwt为二次型性能指标的控制量加权矩阵R；
%  M为控制时域长度；P为预测时域长度，P=Inf时，表示无限的预测和控制长度；
%  tend为仿真结束时间；r为输出设定值或者参考轨迹
%  之后的输入参数为可选参数

%% 进行模型预测控制器设计
plant=model; 
p=6; 
m=2;
ywt=[]; % Q
uwt=[1 1]; % R
r=[1 1]; % 参考轨迹
tend=30;
ulim=[-0.1 -0.1 0.5 0.5 0.1 100]; % 输入下界、上界和变化率曲线的轨迹约束
ylim=[];
[y,u,ym]=cmpc(plant,model,ywt,uwt,m,p,tend,r,ulim,ylim);
plotall(y,u,delt)