matlab运算维数不匹配,matlab编程遇到维数不符合和函数调用出错的问题

winner245 发表于 2013-11-19 08:45

请给出完整的代码

function z

clear

clc

%--------------求稳态值--------------%

us=-10:0.1:10;

ys(1)=0;

for i=1:200;

for k=2:200;

ys(k)=ys(k-1)/(1+ys(k-1)^2)+us(i)^3;

end

YS(i,1)=ys(k);

US(i,1)=us(i);

end   %求出稳态值

figure(1);

plot(US,YS);

title('输入输出稳态非线性特性曲线');

xlabel('输入us');

ylabel('输出ys');

%--------------求稳态模型--------------%

p=polyfit(US,YS,3);%稳态多项式模型系数

KC=polyval(polyder(p),0);%稳态模型在初始点US=0的斜率

%--------------动态模型建立-------------%

y(1)=0.0000;

u(1)=0.1;

for k=2:200;

u(k)=0.1;

y(k)=y(k-1)/(1+y(k-1)^2)+u(k)^3;

Y(k,1)=y(k);

U(k,1)=u(k);

end  %为动态模型辨识提供数据

sys=arx([Y,U],[2,2,1]);

a1=-1.976;

a2=0.9761;

b1=0.0002415;

b2=3.883e-05;

k1=(1+a1+a2)/(b1+b2);

v1=b1*KC*k1;

v2=b2*KC*k1;

%%--------------输出预测-------------%%

P=10;%预测步数

M=5;%控制步数

H1=[0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1];

H2=[0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1];

H3=[0.1,0.2,0.3,0.4,0.5];

A=diag(H1);

B=diag(H3);

N=[0.1 -1.25];%设定反馈增益矩阵0.1z(-1)-1.25

for i=1:P-M    %求C矩阵

for j=1:P

if (j-i==3)&&(i<=P-M)

C(i,j)=H2(M+i)*N(1);

else if (j-i==4)&&(i<=P-M)

C(i,j)=H2(M+i)*N(2);

end

end

end

end           %C矩阵求完

para=[0 0.1 -1.25];%求Q矩阵

r=[0.9761 -1.876 -0.25];

Q=H2(P)*para.'*para+H1(P)*r.'*r; %求得使系统渐进稳定的正定矩阵Q

R=size(r,2);

for i=1:R    %求E矩阵

for j=1:R

if j-i==1

E(i,j)=1;

else if i==R

E(i,j)=-r(j);

end

end

end

end           %E矩阵求完

function F=myfun(D)

F=E'*D*E-D+Q;

end

D0=ones(3,3);

options = optimset('Display','off');

[D,Fval,exitflag] = fsolve(@myfun,D0,options);   %求解D矩阵

for i=1:10

for j=1:10

if 8<=i&&i<=10&&8<=j&&j<=10

DM(i,j)=D(i-7,j-7);

else DM(i,j)=0;

end

end

end  %求矩阵

AM=A+DM;

u(1)=0;%初始值

u(2)=0;%初始值

y(1)=0;%初始值

y(2)=0;%初始值

yds(1:900)=[zeros(1,300),0.3*ones(1,300),0.6*ones(1,300)];%设定值变化

I=1:900;

figure(2);

plot(I,yds)  %画设定值变化曲线

s=0.5;       %设柔化因子

for i=1:P;

SO(i,1)=s^i;

SO1(i,1)=1-s^i;

end%柔化系数

lb=[-3,-3,-3,-3,-3]';

ub=[3,3,3,3,3]';

%%--------------滚动优化------------%%

for k=3:900

YDS=[0 0 0 yds(k)];

pn=p-YDS;

r=roots(pn);

udsf=r(3);

KF=polyval(polyder(p),udsf);%稳态模型在最终稳态点UF的斜率

k2=(KF-KC)/(udsf-0);

g1=b1*k1*k2;

g2=b2*k1*k2;

S(1)=1;

S(2)=a1;

T(1)=v1;

T(2)=a1*v1+v2;

Tg(1)=g1;

Tg(2)=a1*g1+g2;

for n=3:P+1;

S(n)=S(n-1)*a1+S(n-2)*a2;

T(n)=S(n)*v1+S(n-1)*v2;

Tg(n)=S(n)*g1+S(n-1)*g2;

end    %求多步输出预测表达式的系数

for  i=1:P;

G11(i,1)=S(i)*v2;

end  %G11

for j=1:M;

z=0;

for i=1:P;

if i>j-1;

G1(i,j)=T(1+z);

z=z+1;

else

G1(i,j)=0;

end

end

end    %求矢量矩阵G1

for i=1:P;

G22(i,1)=S(i)*g2;

end   %G22

for j=1:M;

z=0;

for i=1:P;

if i>=j-1;

G2(i,j)=Tg(1+z);

z=z+1;

else

G2(i,j)=0;

end

end

end     %求矢量矩阵G2

for i=1:P;

F1(i,1)=S(i+1);

F2(i,1)=a2*S(i);

end  %求矢量矩阵F1和F2

opt=optimset('algorithm','sqp');

U0=ones(M,1);

U=fmincon(@opfun,U0,[],[],[],[],lb,ub,[],opt);

U1(:,k)=U;

u(k)=U1(1,k);

y(k+1)=-a1*y(k)-a2*y(k-1)+v1*u(k)+v2*u(k-1)+g1*u(k)^2+g2*u(k-1)^2;

end

function J=opfun(U)

Y(:,k)=G1*U+G2*U.^2+G11*u(k-1)+G22*u(k-1)^2+F1*y(k)+F2*y(k-1);

Yr(:,k)=SO*y(k)+SO1*yds(k);

J=(Y(:,k)-Yr(:,k)).'*AM*(Y(:,k)-Yr(:,k))+U.'*B*U+norm(C*(Y(:,k)-Yr(:,k)),2)^2;

end

k=1:600;

figure(3);

plot(k,yds(k),'k-',k,y(k),'k:');

end