clear
disp('单变量系统的隐式广义预测控制算法的研究')
disp('广义预测控制算法初始值')
nn=input('时域长度nn=');
n=input('预测长度n=');
m=input('控制长度m=');
t0=input('控制加权系数t0=');
a=input('柔化系数a=');
disp('最小二乘公式初始值')
t1=1; %input('遗忘因子=');
d1=input('n+1阶方阵P的形式:0-对角阵,1-方阵:');
d2=input('n+1阶方阵P的初始值:0-自动赋值1e+5,1-键盘输入:');
if(d1==1)
if(d2==1)
P=input('在方括号【】中,输入(n+1)阶方阵P的值:');
else
p=(1e+5)*ones(n+1);
end
else
if(d2==1)
PP=input('在方括号【】中,输入(n+1)阶对角阵P对角线上的值:');
P=diag(PP);
else
P=(1e+5)*eye(n+1);
end
end
%参数初始值
uuu=0;yyy=0;
uu=zeros(n,1);u=zeros(m,1);
yy=zeros(n,1);y1=zeros(n,1);
Q=zeros(n+1,1);Q(1,1)=1;Q(n+1,1)=1;
%产生周期为100,时间为T,幅值为1的方波信号的给定值
T=300;[yr0,t]=gensig('square',100,T,1);
d3=input('输出曲线是否去掉前100步:0-不,1-去掉:');
nm=length(t);
for ij=2:nm
yr=yr0(ij)+1; %产生周期为100,时间为T,幅值在1和2之间变化的方波信号的给定值
%根据系统模型,计算k时刻的输出值y(k)
y=1.496585*yy(n,1)-0.496585*yy(n-1,1)+0.5*uu(n-1,1);
%(nn=6;n=6;m=2;t0=0.8;a=0.3;t1=1) %example9-1仿真模型
%y=2.001676*yy(n,1)-1.24339*yy(n-1,1)+0.24171*yy(n-2,1)+0.23589*uu(n,1);
%example9-2仿真模型
%产生均匀分布的白噪声
a9=0;
for i=1:1
a9=a9+rand;
end
a8=0.01*(a9-6);
%保存k时刻及以前的n个输出值y(k),y(k-1),...,y(k-n),以供模型运算
for i=1:n-1
yy(i,1)=yy(i+1,1);
end
yy(n,1)=y;
yyy=[yyy;y]; %保存各k时刻nm个输出以便绘图
%根据最小二乘公式,由y(k)计算G阵的各元素值g0,g1,...gn
for i=1:n
X(1,i)=uu(i,1);
end
X(1,n+1)=1;
K=P*X'*inv(t1+X*P*X');
P=(eye(n+1)-K*X)*P/t1;
Q=Q+K*(y-X*Q);
%根据元素值g0,g1,...,gn,求G阵
for j=1:m
for i=n:-1:j
i1=n-i+j;
G(i1,j)=Q(i,1);
end
end
%求nn维y0向量(y1为上一时刻的y0向量)
y0(1:nn-1,1)=y1(2:nn,1);y0(nn,1)=y1(nn,1);
y0=y0+(y-y1(1,1));
for i=1:n
y1(i,1)=y0(i,1)+u(1,1)*Q(n+1-i,1);
end
for i=n+1:nn
y1(i,1)=y1(n,1);
end
%根据y0,求n维f向量f(k+1),...f(k+n)
f(1:n,1)=y0(1:n,1);
%由当前k时刻的输出值y(k)和给定值yr,求k时刻以后的n个参考轨迹w(k+1),...w(k+n)
w=a*y+(1-a)*yr;
for i=2:n
w=[w;a^i*y+(1-a^i)*yr];
end
%计算k时刻及以后的m个控制增量Du(k),,,,Du(k+m)
u=inv(G'*G+t0*eye(m))*G'*(w-f);
%保存k时刻及以后的n个控制增量,以供模型运算
for i=1:n-1
uu(i,1)=uu(i+1,1);
end
uu(n,1)=u(1,1);
uuu=[uuu;u(1,1)]; %保存各k时刻的nm个控制增量以便绘图
%控制量限幅
if(u(1,1)>1)
u(1,1)=1;
end
if(u(1,1)<-1)
u(1,1)=-1;
end
end
%绘制给定值、输出值和控制增量曲线
if(d3==1)
%绘制去掉前100步的给定值、输出值和控制增量曲线
yyy1(1:(T-100),1)=yyy(101:T,1);
uuu1(1:(T-100),1)=uuu(101:T,1);
t1(1:(T-100),1)=t(101:T,1)-100;
yr01(1:(T-100),1)=yr0(101:T,1);
subplot(2,1,1):plot(t1,(yr01+1),t1,yyy1);
axis([0,nm-100,0,2.5]);
xlabel('t');ylabel('yr,y')
subplot(2,1,2):plot(t1,uuu1);
axis([0,nm-100,-1.5,1.5]);
xlabel('t');ylabel('u')
else
%绘制完整的给定值,输出值和控制增量曲线
subplot(2,1,1);plot(t,(yr0+1),t,yyy);
axis([0,nm,0,2.5]);
xlabel('t');ylabel('yr,y')
subplot(2,1,2);plot(t,uuu);
axis([0,nm,-1.5,1.5]);
xlabel('t');ylabel('u')
end
Matlab科研辅导帮
- 粉丝: 3w+
- 资源: 7785
最新资源
- 2_1113重卡电动化创新聚焦日-录音转文字.docx
- “互联网+政务服务”整体建设方案 .docx
- 基于SSH的校园管理系统
- 精选微信小程序源码:面包旅行小程序(旅游类小清新版本)小程序(含源码+源码导入视频教程&文档教程,亲测可用)
- ArcGIS Pro SKD - ADGeoDatabase.daml
- C# winform自定义图片控件.zip,拖拽移动,滚轮缩放
- 基于python的dlib库的人脸识别实现
- ArcGIS Pro SDK - ADCore.daml
- rocketmq的客户端
- 精选微信小程序源码:户外旅游小程序(旅游类)小程序(含源码+源码导入视频教程&文档教程,亲测可用)
资源上传下载、课程学习等过程中有任何疑问或建议,欢迎提出宝贵意见哦~我们会及时处理!
点击此处反馈