clear;
disp('初始值')
nn=input('时域长度p=');n=input('预测长度n=');m=input('控制长度m=');t0=input('控制加权系数λ=');a=input('柔化系数α=');
%最小二乘公式初始化
t1=1;
P=(1e+5)*eye(2*n+1);
%参数初始值
uuu1=0;yyy1=0;uuu2=0;yyy2=0;
uu1=zeros(n,1);u1=zeros(m,1);uu2=zeros(n,1);u2=zeros(m,1);
yy1=zeros(n,1);y11=zeros(n,1);yy2=zeros(n,1);y12=zeros(n,1);
Q1=zeros(2*n+1,1);Q1(1,1)=1;Q1(n+1,1)=1;Q1(2*n+1,1)=1;Q2=Q1;
%产生周期为100,时间为T,幅值为1的方波信号的给定值
T=300;[yr0,t]=gensig('square',100,T,1);
d3=0;d3=input('去掉前100输入1;否则回车');
nm=length(t);
for ij=2:nm
yr1=yr0(ij)+1;
yr2=yr1;
y1=(1+exp(-.7))*yy1(n,1)-exp(-.7)*yy1(n-1,1)+0.3*uu1(n,1)+0.2*uu2(n-1,1);
y2=2.0017*yy2(n,1)-1.2434*yy2(n-1,1)+0.2417*yy2(n-2,1)+0.13*uu1(n-1,1)+0.106*uu2(n,1);
%保存k时刻以前的n个输出值
for i=1:n-1
yy1(i,1)=yy1(i+1,1);yy2(i,1)=yy2(i+1,1);
end
yy1(n,1)=y1;yy2(n,1)=y2;
yyy1=[yyy1;y1];yyy2=[yyy2;y2];
%计算G11,G12,G21,G22各元素g0...
for i=1:n
X(1,i)=uu1(i,1);X(1,i+n)=uu2(i,1);
end
X(1,2*n+1)=1;K=P*X'*inv(t1+X*P*X');
P=(eye(2*n+1)-K*X)*P/t1;
Q1=Q1+K*(y1-X*Q1);
Q2=Q2+K*(y2-X*Q2);
%计算各阵
for j=1:m
for i=n:-1:j
i1=n-i+j;
G11(i1,j)=Q1(i,1);G12(i1,j)=Q1(i+n,1);
G21(i1,j)=Q2(i,1);G22(i1,j)=Q2(i+n,1);
end
end
%for i=1:n,G12(i,2)=0;G21(i,2)=0;end
%求nn维y01,y02(y11和y12为上一时刻的y01,y02)
e1=y1-y11(1,1);e2=y2-y12(1,1);
y01(1:nn-1,1)=y11(2:nn,1);y01(nn,1)=y11(nn,1);y01=y01+e1;
y02(1:nn-1,1)=y12(2:nn,1);y02(nn,1)=y12(nn,1);y02=y02+e2;
for i=1:n
y11(i,1)=y01(i,1)+u1(1,1)*Q1(n+1-i,1)+u2(1,1)*Q1(2*n+1-i,1);
y12(i,1)=y02(i,1)+u2(1,1)*Q2(n+1-i,1)+u1(1,1)*Q2(2*n+1-i,1);
end
for i=n+1:nn
y11(i,1)=y11(n,1);y12(i,1)=y12(n,1);
end
%求f向量
f1(1:n,1)=y01(1:n,1);f2(1:n,1)=y02(1:n,1);
%保存k时刻及以后的n个控制增量
for i=1:n-1
uu1(i,1)=uu1(i+1,1);uu2(i,1)=uu2(i+1,1);
end
uu1(n,1)=u1(1,1);uu2(n,1)=u2(1,1);
uuu1=[uuu1,u1(1,1)];uuu2=[uuu2,u2(1,1)];
%求K 时刻以后的n个参考轨迹
w1=a*y1+(1-a)*yr1;w2=a*y2+(1-a)*yr2;
for i=2:n
w1=[w1;a^i*y1+(1-a^i)*yr1];w2=[w2;a^i*y2+(1-a^i)*yr2];
end
%计算k时刻以后m个控制增量
u1=inv(G11'*G11+t0*eye(m))*G11'*(w1-G12*u2-f1);
u2=inv(G22'*G22+t0*eye(m))*G22'*(w2-G21*u1-f2);
%控制量限幅
if(u1(1,1)>1),u1(1,1)=1;end
if(u1(1,1)<-1),u1(1,1)=-1;end
if(u2(1,1)>1),u2(1,1)=1;end
if(u2(1,1)<-1),u2(1,1)=-1;end
end
%绘制曲线
if(d3==1)
yyyy1(1:(T-100),1)=yyy1(101:T,1);uuuu1(1:(T-100),1)=uuu1(101:T,1);
yyyy2(1:(T-100),1)=yyy2(101:T,1);uuuu2(1:(T-100),1)=uuu2(101:T,1);
t1(1:(T-100),1)=t(101:T,1)-100;yr01(1:(T-100),1)=yr0(101:T,1);
subplot(2,2,1):plot(t1,(yr01+1),t1,yyyy1);axis([0,nm-100,0,2.5]);xlabel('t');ylabel('yr1,y1');
subplot(2,2,3):plot(t1,uuuu1);axis([0,nm-100,-1.5,1.5]);xlabel('t');ylabel('u1');
subplot(2,2,2):plot(t1,(yr01+1),t1,yyyy2);axis([0,nm-100,0,2.5]);xlabel('t');ylabel('yr2,y2');
subplot(2,2,4):plot(t1,uuuu2);axis([0,nm-100,-1.5,1.5]);xlabel('t');ylabel('u2');
else
subplot(2,2,1):plot(t,(yr0+1),t,yyy1);axis([0,nm,0,2.5]);xlabel('t');ylabel('yr1,y1');
subplot(2,2,3):plot(t,uuu1);axis([0,nm,-1.5,1.5]);xlabel('t');ylabel('u1');
subplot(2,2,2):plot(t,(yr0+1),t,yyy2);axis([0,nm,0,2.5]);xlabel('t');ylabel('yr2,y2');
subplot(2,2,4):plot(t,uuu1);axis([0,nm,-1.5,1.5]);xlabel('t');ylabel('u2');
end
广义预测控制
需积分: 18 19 浏览量
2018-03-14
16:04:05
上传
评论 1
收藏 2KB RAR 举报
weixin_41839959
- 粉丝: 0
- 资源: 1
最新资源
- AIS2024 valid
- 最入门的爬虫代码 python.docx
- 爬虫零基础入门-爬取天气预报.pdf
- 最通俗易懂的 MongoDB 非结构化文档存储数据库教程.zip
- 以mongodb为数据库的订单物流小项目.zip
- 腾讯云-mongodb数据库, 项目部署.zip
- 腾讯 APIJSON 的 MongoDB 数据库插件.zip
- 理解非关系型数据库和关系型数据库的区别.zip
- 操作简单的Mongodb网页web管理工具,基于Spring Boot2.0支持mongodb集群.zip
- tms-mongodb-web,提供访问mongodb数据的REST API和可灵活扩展的mongodb web 客户端.zip
资源上传下载、课程学习等过程中有任何疑问或建议,欢迎提出宝贵意见哦~我们会及时处理!
点击此处反馈