matlab 可以求稳态解吗,基于Matlab/Simulink的发动机稳态求解

H=0;

Ma=0;

Wf=2.5;

n=3;

k1=0.0017;

k2=0.00173;

k3=0.002;

k4=0.00882;

k5=0.00345;

k6=0.00516;

k7=10;

k8=4;

Jacobian=zeros(8);

A=[PFc,PFd,Pbst,Phpc,Phpt,Plpt,N2,N1];

ZB=zeros(1,n);

%ZB=sqrt(f1^2+f2^2+f3^2+f4^2+f5^2+f6^2+f7^2+f8^2);

for i=1:1:n

X_guess=A;

[t,x,y]=sim('allengine',0);

f1=y(:,18);

f2=y(:,19);

f3=y(:,20);

f4=y(:,21);

f5=y(:,22);

f6=y(:,23);

f7=y(:,24);

f8=y(:,25);

ZB(1,i)=(f1^2+f2^2+f3^2+f4^2+f5^2+f6^2+f7^2+f8^2)^(0.5);

%以下针对第一个变量进行扰动

B1=[A(1)+k1,A(2),A(3),A(4),A(5),A(6),A(7),A(8)];

X_guess=B1;

[t,x,y]=sim('allengine',0);

f11=y(:,18);

Jacobian(1,1)=(f11-f1)/(k1);

f21=y(:,19);

Jacobian(2,1)=(f21-f2)/(k1);

f31=y(:,20);

Jacobian(3,1)=(f31-f3)/(k1);

f41=y(:,21);

Jacobian(4,1)=(f41-f4)/(k1);

f51=y(:,22);

Jacobian(5,1)=(f51-f5)/(k1);

f61=y(:,23);

Jacobian(6,1)=(f61-f6)/(k1);

f71=y(:,24);

Jacobian(7,1)=(f71-f7)/(k1);

f81=y(:,25);

Jacobian(8,1)=(f81-f8)/(k1);

%以下是针对第二个变量

B2=[A(1),A(2)+k2,A(3),A(4),A(5),A(6),A(7),A(8)];

X_guess=B2;

[t,x,y]=sim('allengine',0);%%运行simulink%%

f12=y(:,18);

Jacobian(1,2)=(f12-f1)/(k2);

f22=y(:,19);

Jacobian(2,2)=(f22-f2)/(k2);

f32=y(:,20);

Jacobian(3,2)=(f32-f3)/(k2);

f42=y(:,21);

Jacobian(4,2)=(f42-f4)/(k2);

f52=y(:,22);

Jacobian(5,2)=(f52-f5)/(k2);

f62=y(:,23);

Jacobian(6,2)=(f62-f6)/(k2);

f72=y(:,24);

Jacobian(7,2)=(f72-f7)/(k2);

f82=y(:,25);

Jacobian(8,2)=(f82-f8)/(k2);

%以下是针对第三个变量

B3=[A(1),A(2),A(3)+k3,A(4),A(5),A(6),A(7),A(8)];

X_guess=B3;

[t,x,y]=sim('allengine',0);%%运行simulink%%

f13=y(:,18);

Jacobian(1,3)=(f13-f1)/(k3);

f23=y(:,19);

Jacobian(2,3)=(f23-f2)/(k3);

f33=y(:,20);

Jacobian(3,3)=(f33-f3)/(k3);

f43=y(:,21);

Jacobian(4,3)=(f43-f4)/(k3);

f53=y(:,22);

Jacobian(5,3)=(f53-f5)/(k3);

f63=y(:,23);

Jacobian(6,3)=(f63-f6)/(k3);

f73=y(:,24);

Jacobian(7,3)=(f73-f7)/(k3);

f83=y(:,25);

Jacobian(8,3)=(f83-f8)/(k3);

%以下是针对第四个变量

B4=[A(1),A(2),A(3),A(4)+k4,A(5),A(6),A(7),A(8)];

X_guess=B4;

[t,x,y]=sim('allengine',0);%%运行simulink%%

f14=y(:,18);

Jacobian(1,4)=(f14-f1)/(k4);

f24=y(:,19);

Jacobian(2,4)=(f24-f2)/(k4);

f34=y(:,20);

Jacobian(3,4)=(f34-f3)/(k4);

f44=y(:,21);

Jacobian(4,4)=(f44-f4)/(k4);

f54=y(:,22);

Jacobian(5,4)=(f54-f5)/(k4);

f64=y(:,23);

Jacobian(6,4)=(f64-f6)/(k4);

f74=y(:,24);

Jacobian(7,4)=(f74-f7)/(k4);

f84=y(:,25);

Jacobian(8,4)=(f84-f8)/(k4);

%以下是针对第五个变量

B5=[A(1),A(2),A(3),A(4),A(5)+k5,A(6),A(7),A(8)];

X_guess=B5;

[t,x,y]=sim('allengine',0);%%运行simulink%%

f15=y(:,18);

Jacobian(1,5)=(f15-f1)/(k5);

f25=y(:,19);

Jacobian(2,5)=(f25-f2)/(k5);

f35=y(:,20);

Jacobian(3,5)=(f35-f3)/(k5);

f45=y(:,21);

Jacobian(4,5)=(f45-f4)/(k5);

f55=y(:,22);

Jacobian(5,5)=(f55-f5)/(k5);

f65=y(:,23);

Jacobian(6,5)=(f65-f6)/(k5);

f75=y(:,24);

Jacobian(7,5)=(f75-f7)/(k5);

f85=y(:,25);

Jacobian(8,5)=(f85-f8)/(k5);

%以下是针对第六个变量

B6=[A(1),A(2),A(3),A(4),A(5),A(6)+k6,A(7),A(8)];

X_guess=B6;

[t,x,y]=sim('allengine',0);%%运行simulink%%

f16=y(:,18);

Jacobian(1,6)=(f16-f1)/(k6);

f26=y(:,19);

Jacobian(2,6)=(f26-f2)/(k6);

f36=y(:,20);

Jacobian(3,6)=(f36-f3)/(k6);

f46=y(:,21);

Jacobian(4,6)=(f46-f4)/(k6);

f56=y(:,22);

Jacobian(5,6)=(f56-f5)/(k6);

f66=y(:,23);

Jacobian(6,6)=(f66-f6)/(k6);

f76=y(:,24);

Jacobian(7,6)=(f76-f7)/(k6);

f86=y(:,25);

Jacobian(8,6)=(f86-f8)/(k6);

%以下是针对第七个变量

B7=[A(1),A(2),A(3),A(4),A(5),A(6),A(7)+k7,A(8)];

X_guess=B7;

[t,x,y]=sim('allengine',0);%%运行simulink%%

f17=y(:,18);

Jacobian(1,7)=(f17-f1)/(k7);

f27=y(:,19);

Jacobian(2,7)=(f27-f2)/(k7);

f37=y(:,20);

Jacobian(3,7)=(f37-f3)/(k7);

f47=y(:,21);

Jacobian(4,7)=(f47-f4)/(k7);

f57=y(:,22);

Jacobian(5,7)=(f57-f5)/(k7);

f67=y(:,23);

Jacobian(6,7)=(f67-f6)/(k7);

f77=y(:,24);

Jacobian(7,7)=(f77-f7)/(k7);

f87=y(:,25);

Jacobian(8,7)=(f87-f8)/(k7);

%以下是针对第八个变量

B8=[A(1),A(2),A(3),A(4),A(5),A(6),A(7),A(8)+k8];

X_guess=B8;

[t,x,y]=sim('allengine',0);%%运行simulink%%

f18=y(:,18);

Jacobian(1,8)=(f18-f1)/(k8);

f28=y(:,19);

Jacobian(2,8)=(f28-f2)/(k8);

f38=y(:,20);

Jacobian(3,8)=(f38-f3)/(k8);

f48=y(:,21);

Jacobian(4,8)=(f48-f4)/(k8);

f58=y(:,22);

Jacobian(5,8)=(f58-f5)/(k8);

f68=y(:,23);

Jacobian(6,8)=(f68-f6)/(k8);

f78=y(:,24);

Jacobian(7,8)=(f78-f7)/(k8);

f88=y(:,25);

Jacobian(8,8)=(f88-f8)/(k8);

M=pinv(Jacobian);

N=M*[f1;f2;f3;f4;f5;f6;f7;f8];  %逆矩阵乘以原误差值矩阵

A(1)=A(1)-N(1,1);  %最新的n_L等参数

A(2)=A(2)-N(2,1);

A(3)=A(3)-N(3,1);

A(4)=A(4)-N(4,1);

A(5)=A(5)-N(5,1);

A(6)=A(6)-N(6,1);

A(7)=A(7)-N(7,1);

A(8)=A(8)-N(8,1);

end

B=A;

%C=zeros(1,n);

%for i=1:n

%  C(i)=i;

%end

%plot(C,ZB,'.');xlabel('迭代次数');ylabel('方程组解的误差');

%title('设计点共同工作图');

%set(gca, 'fontsize', 12);  %更改坐标轴字体大小

%}