预测控制作业程序调试
8-1：考虑一个单输入单输出的对象，其传递函数为：
2
1
()36
s G s s s
+=
++
解：
采用下面的MATLAB 程序对该对象进行脉冲响应模型辨识,脉冲响应模型预测输出与预测误差曲线：
% 将多项式的传递函数模型转换为MPC 传递函数模型 num=[1 1];den=[1 3 6]; h=tf(num,den); % 获得脉冲信号x
[u,t]=gensig('pulse',2,10,0.1); x=u;
% 求解LTI 对象的单位脉冲响应y t=0:0.1:10;
[y,x,t]=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); 结果如图所示：
8-2：考虑一个双输入单输出的对象，其传递函数矩阵为：
14155.72 1.52601251
()s s
e e s s G s --++??
=?
?
采样时间为7秒 解：
采用下面的MATLAB 程序对该对象进行脉冲响应模型辨识,脉冲响应模型预测输出与预测误差曲线： 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 传递函数模型转换为MPC 状态空间模型 mod=tfd2mod(7,1,g1,g2);
%将MPC 状态空间模型转换为通用状态空间模型 [A,B,C,D]=mod2ss(mod);
%将通用状态空间模型转换LTI 对象的 状态空间模型 sys=ss(A,B,C,D); h=tf(sys);
% 获得脉冲信号x
[u,t]=gensig('pulse',8,10,0.1);
x=[u u];
% 求解LTI 对象的单位脉冲响应y t=0:0.1:10;
[y,x1,t1]=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);
% 基于多变量最小二乘法的脉冲响应模型辨识 ninput=2; plotopt=2;
[theta,yres]=mlr(xreg,yreg,ninput,plotopt); 结果如图所示：
8-3：考虑一个单输入单输出的对象，其传递函数为： 2
1
()36
s G s s s
+=
++
解：
采用下面的MATLAB 程序对该对象进行脉冲响应模型辨识,脉冲响应模型预测输出与预测误差曲线
% 将多项式的传递函数模型转换为MPC 传递函数模型 num=[1 1];den=[1 3 6]; h=tf(num,den); % 获得脉冲信号x
[u,t]=gensig('pulse',2,10,0.1); x=[u];
% 求解LTI 对象的单位脉冲响应y t=0:0.1:10;
[y,x1,t1]=lsim(h,x,t);
% 生成用于线性回归计算的输入/输出数据矩阵 n=30;
[xreg,yreg] = wrtreg(Sx,y,n);
% 基于部分最小二乘法的脉冲响应模型辨识 ninput=2; lv=10; plotopt=2;
[theta,yres]=mlr(xreg,yreg,ninput,plotopt); 结果如图所示：
8-4：考虑一个单输入单输出的对象，其传递函数为： 2
1
()36
s G s s s
+=
++
解：
采用下面的MATLAB 程序对该对象进行脉冲响应模型辨识,脉冲响应模型预测输出与预测误差曲线： num=[1 1];den=[1 3 6]; h=tf(num,den); % 获得脉冲信号x
[u,t]=gensig('pulse',2,10,0.1); x=[u];
% 求解LTI 对象的单位脉冲响应y t=0:0.1:10;
[y,x1,t1]=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)； 结果如图所示：
8-12：考虑如下的双输入双输出纯时延对象，其传递函数距阵为：
73312.8 6.616.7110.9118.919.4211
14.41
()s s s
s e e s s e e s s G s ----++--++?
?= ?
???
解：
MATLAB 程序如下：
%将传递函数模型转换为阶跃响应模型
g11=poly2tfd(12.8,[16.7 1],0,1); g12 = poly2tfd(6.6,[10.9 1],0,7); g21 = poly2tfd(-18.9,[21.0 1],0,3); g22 = poly2tfd(-19.4,[14.4 1],0,3); delt=3; ny=2; tfinal=90;
model = tfd2step(tfinal, delt, ny, g11, g12,g21,g22); %进行模型预测控制器设计 plant=model;
%预测时域长度为6 p=6;m=2;
ywt=[];uwt=[1 1];
%设置输入约束和参考轨迹等控制器参数 r=[1 1];
tend=30;%仿真时间为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)
闭环系统的输出和控制量变化曲线如图所示：
8-15：设系统的传递函数距阵为：
73312.8 6.616.7110.9118.919.4211
14.41
()s s s
s e e s s e e s s G s ----++--++?
?= ?
???
解：
MATLAB 程序如下：
g11=poly2tfd(12.8,[16.7 1],0,1); g12 = poly2tfd(6.6,[10.9 1],0,7); g21 = poly2tfd(-18.9,[21.0 1],0,3); g22 = poly2tfd(-19.4,[14.4 1],0,3); delt=3; ny=2;
imod = tfd2mod(delt, ny, g11, g12,g21,g22); pmod=imod; p=6;m=2;
ywt=[];uwt=[1 1]; tend=30; r=[0 1];
ulim=[-inf -0.15 inf inf 0.1 100]; ylim=[];
[y,u,ym]=scmpc(pmod,imod,ywt,uwt,m,p,tend,r,ulim,ylim); plotall(y,u,delt)
闭环系统的输出和控制量变化曲线如图所示：
在添加对输出变量的约束后，再利用以下程序对系统进行模型预测控制器的设计，得到闭环控制系统输出响应和控制量变化曲线如图： ulim=[-inf -0.15 inf inf 0.1 100]; ylim=[0 0 inf inf];
[y,u,ym]=scmpc(pmod,imod,ywt,uwt,m,p,tend,r,ulim,ylim);
plotall(y,u,delt)
8-16：考虑有如下传递函数距阵的多变量系统的状态空间模型预测控制器设计问题
73312.8 6.616.7110.9118.919.4211
14.41
()s s s
s e e s s e e s s G s ----++--++?
?= ? ???
解：
MATLAB 程序如下：
%在进行模型预测控制器设计之前，首先将系统模型转换为状态空间形式
T=2;
g11=poly2tfd(12.8,[16.7 1],0,1);
g12 = poly2tfd(6.6,[10.9 1],0,7);
g21 = poly2tfd(-18.9,[21.0 1],0,3);
g22 = poly2tfd(-19.4,[14.4 1],0,3);
umod = tfd2mod(T,2, g11, g12,g21,g22);
%定义扰动模型
g13=poly2tfd(3.8,[14.9 1],0,8)
g23=poly2tfd(4.9,[13.2 1],0,3)
dmod=tfd2mod(T,2,g13,g23);
%建立叠加了扰动的混合系统模型
pmod=addumd(umod,dmod);
%考虑精确建模的情况
imod=pmod;
ywt=[];uwt=[];
%预测时域和控制时域均为5
P=5;
M=P;
Ks=smpccon(imod,ywt,uwt,M,P);
tend=30;
r=[0 1];
[y,u]=smpcsim(pmod,imod,Ks,tend,r);
plotall(y,u,T)
得到闭环控制系统输出响应和控制量变化曲线如图：
增加预测时域长度，同时减少控制时域长度后，再利用以下程序可得闭环系统的输出和控制量变化曲线如图所示：
P=10;
M=3;
Ks=smpccon(imod,ywt,uwt,M,P);
[y,u]=smpcsim(pmod,imod,Ks,tend,r);
plotall(y,u,T)
进一步改变控制时域长度，采用控制量分块的形式，即再利用以下程序可得到的闭环系统的输出和控制量变化曲线
M=[2 3 4];
Ks=smpccon(imod,ywt,uwt,M,P);
[y,u]=smpcsim(pmod,imod,Ks,tend,r);
plotall(y,u,T)
增加输入控制量的加权距阵系数，即再利用以下程序可得模型预测闭环控制系统输出控制量变化曲线
uwt=[1 1];
P=5;
M=P;
Ks=smpccon(imod,ywt,uwt,M,P);
tend=30;
r=[1 0];
[y,u]=smpcsim(pmod,imod,Ks,tend,r);
plotall(y,u,T)
再利用以下程序将输出设定值均设为0，绘制闭环系统的输出和控制量变化曲线，如图所示：
ulim=[];Kest=[];R=[];
z=[];
v=[];
w=[1];
[y,u]=smpcsim(pmod,imod,Ks,tend,r,ulim,Kest,z,v,w);
plotall(y,u,T)
利用以下程序采用估计器进一步改善系统性能，对应的系统输出和控制量变化曲
线如图所示
[Kest,newmod]=smpcest(imod,[15,15],[3 3]); Ks=smpccon(newmod,ywt,uwt,M,P);
[y,u]=smpcsim(pmod,newmod,Ks,tend,r,ulim,Kest,z,v,w);
plotall(y,u,T)
8-17：考虑如下的多变量系统
73312.8 6.616.7110.9118.919.4211
14.41
()s s s
s e e s s e e s s G s ----++--++?
?= ?
???
解：
MATLAB 程序如下：
g11=poly2tfd(12.8,[16.7 0],0,1); g12 = poly2tfd(6.6,[10.9 1],0,7); g21 = poly2tfd(-18.9,[21.0 1],0,3); g22 = poly2tfd(-19.4,[14.4 1],0,3); delt=1; ny=2;
imod = tfd2mod(delt,ny, g11, g12,g21,g22); gw1=poly2tfd(3.8,[14.9 1],0,8) gw2=poly2tfd(4.9,[13.2 1],0,3)
pmod=addumd(imod,tfd2mod(delt,ny,gw1,gw2)); %设计模型预测控制器 P=6;M=2;
ywt=[];uwt=[1 1];
Ks=smpccon(imod,ywt,uwt,M,P);
r=[];
ulim=[];
z=[];
v=[];
w=[1];
wu=[];
tend=30;
[y3,u3]=smpcsim(pmod,imod,Ks,tend,r,ulim,[],z,v,w,wu);
%设计状态估计器
Kest1=smpcest(pmod,1,0.001*eye(ny));
Ks1=smpccon(pmod,ywt,uwt,M,P);
[y1,u1]=smpcsim(pmod,pmod,Ks1,tend,r,ulim,Kest1,z,v,w,wu); plotall(y1,u1,delt)
采用了状态估计器的闭环系统输出和控制量变化曲线如图所示：
下面进行简化的状态估计器设计
tau=[10 10];
signoise=[3 3];
[Kest2,newmod]=smpcest(imod,tau,signoise);
Ks2=smpccon(newmod,ywt,uwt,M,P);
[y2,u2]=smpcsim(pmod,pmod,Ks1,tend,r,ulim,Kest1,z,v,w,wu); plotall(y2,u2,delt)
采用简化的状态估计器的闭环系统输出和控制量曲线如图所示：
专业：控制理论与控制工程 学号：062030040 姓名：陈孝凯
程序调试
8-1：考虑一个单输入单输出的对象，其传递函数为：
2
1
()36
s G s s s
+=
++
解：
采用下面的MATLAB 程序对该对象进行脉冲响应模型辨识,脉冲响应模型预测输出与预测误差曲线：
% 将多项式的传递函数模型转换为MPC 传递函数模型 num=[1 1];den=[1 3 6]; h=tf(num,den); % 获得脉冲信号x
[u,t]=gensig('pulse',2,10,0.1); x=u;
% 求解LTI 对象的单位脉冲响应y t=0:0.1:10;
[y,x,t]=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); 结果如图所示：
8-2：考虑一个双输入单输出的对象，其传递函数矩阵为：
14155.72 1.52601251
()s s
e e s s G s --++??
=?
?
采样时间为7秒 解：
采用下面的MATLAB 程序对该对象进行脉冲响应模型辨识,脉冲响应模型预测输出与预测误差曲线： 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 传递函数模型转换为MPC 状态空间模型 mod=tfd2mod(7,1,g1,g2);
%将MPC 状态空间模型转换为通用状态空间模型 [A,B,C,D]=mod2ss(mod);
%将通用状态空间模型转换LTI 对象的 状态空间模型 sys=ss(A,B,C,D); h=tf(sys);
% 获得脉冲信号x
[u,t]=gensig('pulse',8,10,0.1);
x=[u u];
% 求解LTI 对象的单位脉冲响应y t=0:0.1:10;
[y,x1,t1]=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);
% 基于多变量最小二乘法的脉冲响应模型辨识 ninput=2; plotopt=2;
[theta,yres]=mlr(xreg,yreg,ninput,plotopt); 结果如图所示：
8-3：考虑一个单输入单输出的对象，其传递函数为： 2
1
()36
s G s s s
+=
++
解：
采用下面的MATLAB 程序对该对象进行脉冲响应模型辨识,脉冲响应模型预测输出与预测误差曲线
% 将多项式的传递函数模型转换为MPC 传递函数模型 num=[1 1];den=[1 3 6]; h=tf(num,den); % 获得脉冲信号x
[u,t]=gensig('pulse',2,10,0.1);
x=[u];
% 求解LTI 对象的单位脉冲响应y t=0:0.1:10;
[y,x1,t1]=lsim(h,x,t);
% 生成用于线性回归计算的输入/输出数据矩阵 n=30;
[xreg,yreg] = wrtreg(Sx,y,n);
% 基于部分最小二乘法的脉冲响应模型辨识 ninput=2; lv=10; plotopt=2;
[theta,yres]=mlr(xreg,yreg,ninput,plotopt); 结果如图所示：
8-4：考虑一个单输入单输出的对象，其传递函数为： 2
1
()36
s G s s s
+=
++
解：
采用下面的MATLAB 程序对该对象进行脉冲响应模型辨识,脉冲响应模型预测输出与预测误差曲线： num=[1 1];den=[1 3 6]; h=tf(num,den); % 获得脉冲信号x
[u,t]=gensig('pulse',2,10,0.1); x=[u];
% 求解LTI 对象的单位脉冲响应y t=0:0.1:10;
[y,x1,t1]=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)； 结果如图所示：
8-12：考虑如下的双输入双输出纯时延对象，其传递函数距阵为：
73312.8 6.616.7110.9118.919.4211
14.41
()s s s
s e e s s e e s s G s ----++--++?
?= ?
???
解：
MATLAB 程序如下：
%将传递函数模型转换为阶跃响应模型
g11=poly2tfd(12.8,[16.7 1],0,1); g12 = poly2tfd(6.6,[10.9 1],0,7); g21 = poly2tfd(-18.9,[21.0 1],0,3);
delt=3; ny=2; tfinal=90;
model = tfd2step(tfinal, delt, ny, g11, g12,g21,g22); %进行模型预测控制器设计 plant=model;
%预测时域长度为6 p=6;m=2;
ywt=[];uwt=[1 1];
%设置输入约束和参考轨迹等控制器参数 r=[1 1];
tend=30;%仿真时间为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)
闭环系统的输出和控制量变化曲线如图所示：
8-15：设系统的传递函数距阵为：
73312.8 6.616.7110.9118.919.4211
14.41
()s s s
s e e s s e e s s G s ----++--++?
?= ?
???
解：
MATLAB 程序如下：
g11=poly2tfd(12.8,[16.7 1],0,1); g12 = poly2tfd(6.6,[10.9 1],0,7); g21 = poly2tfd(-18.9,[21.0 1],0,3);
delt=3;
ny=2;
imod = tfd2mod(delt, ny, g11, g12,g21,g22);
pmod=imod;
p=6;m=2;
ywt=[];uwt=[1 1];
tend=30;
r=[0 1];
ulim=[-inf -0.15 inf inf 0.1 100];
ylim=[];
[y,u,ym]=scmpc(pmod,imod,ywt,uwt,m,p,tend,r,ulim,ylim);
plotall(y,u,delt)
闭环系统的输出和控制量变化曲线如图所示：
在添加对输出变量的约束后，再利用以下程序对系统进行模型预测控制器的设计，得到闭环控制系统输出响应和控制量变化曲线如图：
ulim=[-inf -0.15 inf inf 0.1 100];
ylim=[0 0 inf inf];
[y,u,ym]=scmpc(pmod,imod,ywt,uwt,m,p,tend,r,ulim,ylim);
plotall(y,u,delt)
8-16：考虑有如下传递函数距阵的多变量系统的状态空间模型预测控制器设计问题
73312.8 6.616.7110.9118.919.4211
14.41
()s s s
s e e s s e e s s G s ----++--++?
?= ? ???
解：
MATLAB 程序如下：
%在进行模型预测控制器设计之前，首先将系统模型转换为状态空间形式 T=2;
g11=poly2tfd(12.8,[16.7 1],0,1); g12 = poly2tfd(6.6,[10.9 1],0,7); g21 = poly2tfd(-18.9,[21.0 1],0,3); g22 = poly2tfd(-19.4,[14.4 1],0,3);
umod = tfd2mod(T,2, g11, g12,g21,g22); %定义扰动模型
g13=poly2tfd(3.8,[14.9 1],0,8) g23=poly2tfd(4.9,[13.2 1],0,3) dmod=tfd2mod(T,2,g13,g23);
%建立叠加了扰动的混合系统模型 pmod=addumd(umod,dmod); %考虑精确建模的情况 imod=pmod; ywt=[];uwt=[];
%预测时域和控制时域均为5 P=5; M=P;
Ks=smpccon(imod,ywt,uwt,M,P); tend=30;