基于状态模拟方法分析高压油管压力控制问题---源代码（matlab）

应一些读者要求，这里加上了源代码，奈何本人最近忙于数学建模培训，所以排版可能不太友好，但是里面每一个都写了序号，你们可以看一下

1.m1.m
clear;clc;
load data1
y=E;
x=[ones(size(y)) p.^2 p];
[b,bint,r,rint,stats]=regress(y,x)
P=b(1)+b(2)*(p.^2)+b(3)*p;
figure;
plot(p,E,'r--',p,P,'b');
title("tan xing mo liang he ya qiang");
xlabel("ya qiang");
ylabel("tan xing mo liang");
legend('shi ji zhi','hui gui qu xian');

2.m2.m
clear;clc;
rho=0.85:0.005:0.90;
P=-53.23+227.0*tan(1.659+6.558*log(rho));
figure;
plot(rho,P);
title('???????');
xlabel('??(mg/mm^3)');
ylabel('??(MPa)');

3.m3.m
clear;clc;
load data1;
x=100:0.1:150;
y=1./(b(2).*x.^2+b(3).*x+b(1));
rho0=0.85;
rho=exp(trapz(x,y)+log(rho0))
x=100:0.1:160;
y=1./(b(2).*x.^2+b(3).*x+b(1));
rho0=0.85;
rho=exp(trapz(x,y)+log(rho0))
4.m4.m
clear;clc;
options=struct();
options.yfun=@yfun_1;
options.target=@target_1;
options.isplot=0;
options.interval=0.1;
options.duration=2000;
[fm,mx]=stepsolve(@targetfun,0.26,0.31,0.0005,0.02,5,options);
options.isplot=1;
targetfun(mx,options);
ylabel('??(MPa)');
xlabel('??(ms)');
title('???????????');
5.m5.m
clear;clc;
options=struct();
options.yfun=@yfun_1;
options.target=@target_2;
options.isplot=0;
options.interval=1;
options.duration=10000;
[fm,mx]=stepsolve(@targetfun,0.6,1,0.0005,0.05,5,options);
options.isplot=1;
targetfun(mx,options);
ylabel('??(MPa)');
xlabel('??(ms)');
title('???????????');
6.m6.m
clear;clc;
options=struct();
options.yfun=@yfun_2;
options.target=@target_1;
options.isplot=0;
options.interval=1;
options.duration=10000;
[fm,mx]=stepsolve(@targetfun,0.019,0.029,0.000001,0.001,3,options);
options.isplot=1;
targetfun(mx,options);
ylabel('??(MPa)');
xlabel('??(ms)');
title('???????????');
7.m7.m
clear;clc;
options=struct();
options.yfun=@yfun_3;
options.target=@target_1;
options.isplot=0;
options.interval=1;
options.duration=10000;
[fm,mx]=stepsolve(@targetfun,0.039,0.059,0.00001,0.002,3,options);
options.isplot=1;
targetfun(mx,options);
ylabel('高压油管压力(MPa)');
xlabel('运行时间(ms)');
title('两个喷油口工作时高压油管内部压力');
8.m8.m
clear;clc;
options=struct();
options.yfun=@yfun_4;
options.target=@target_1;
options.isplot=0;
options.interval=1;
options.duration=10000;
[fm,mx]=stepsolve(@targetfun,97,103,0.1,1,3,options);
options.isplot=1;
targetfun(mx,options);
ylabel('高压油管压力(MPa)');
xlabel('运行时间(ms)');
title('增加减压阀后高压油管内部压力');
9.Pfun.m
function P=Pfun(x)
if (x<0.7)
P=0;
else
P=(-53.23+227.0*tan(1.659+6.558*log(x)));
end
end
10.Q_in_1.m
function Q=Q_in_1(tA,t,P_in)
%tA????
%t????
%P_in
T=tA+10;
if mod(t,T)<tA
Q=0.85*(pi/4)*(1.4)^2*sqrt(2*(160-P_in)/0.871);
else
Q=0;
end
end
11.Q_in_2.m
function Q=Q_in_2(w,t,P_in)
%tA????
%t????
%P_in
T=tA+10;
if mod(t,T)<tA
Q=0.85*(pi/4)*(1.4)^2*sqrt(2*(160-P_in)/0.871);
else
Q=0;
end
end
12.Q_out_1.m
function Q=Q_out_1(t,P)
r=mod(t,100);
if r<=0.2
Q=100*r;
elseif r<=2.2
Q=20;
elseif r<=2.4
Q=-100*r+240;
else
Q=0;
end
end
13.Q_out_2.m
function Q=Q_out_2(t,p)
%t:时间
%p:密度
r=mod(t,100);
b1=[16.3667705292548,-2.27403710032532,0.0443842537666548];
b2=[16.4024620318978,-78.0705445458018,92.8628884740531];
if r<=0.44
h=polyval(b1,r);
elseif r<=2
h=2;
elseif r<=2.45
h=polyval(b2,r);
else
h=0;
end
dk=1.4;
dz=2.5;
alpha=9/180*pi;
A=min([pi*dk^2/4,pi/4*(4*dz*h*tan(alpha)+2*h^2*tan(alpha)^2)]);
Q=0.85*A*(2*Pfun(p)/p)^0.5;
end
14.Q_out_3.m
function Q= Q_out_3(p,rho)
Q=0.85*0.7^2*pi*sqrt(p/rho);
end
15.stepsolve.m
function [mf,M]=stepsolve(fun,L,R,tol,h0,d,options)
if nargin<=6
options=struct();
end
optfun=@(x)fun(x,options);
M=inf;
fprintf("Initiating...\n");
mf=optfun(M);
while R-L>tol
fprintf("Calculating.. From %.5f to %.5f step %.5f\n",L,R,h0);
t_A=L:h0:R;
ff=arrayfun(optfun,t_A);
mf=min(ff);
M=t_A(mf==ff);
L=M-(M-L)/d;
R=M+(R-M)/d;
h0=h0/d;
fprintf("Current: Min: %.5f At %.5f\n",mf,M);
end
end
16.target_1.m
function [f]=target_1(TT,PP)
f=norm(PP-100,2)/(max(TT)-min(TT));
end
17.target_2.m
function [f]=target_2(TT,PP)
if max(PP)>150
f=sum(diff(TT).*PP(1:end-1));
else
f=1e12;
end
end
18.targetfun.m
function f=targetfun(w,options)
if nargin==1
options=struct();
end
interval=initfun(options,'interval',1);
duration=initfun(options,'duration',20000);
isplot=initfun(options,'isplot',1);
target=options.target;
yfunopt=options.yfun;
m=0;
global flag;
flag=0;
yfun=@(t,y,m)yfunopt(t,y,m,w,interval);
y=0.850;
steps=0:interval:duration;
TT=zeros(size(steps));
YY=zeros(size(steps));
count=0;
for t=steps
count=count+1;
[y,m]=yfun(t,y,m);
TT(count)=t;
YY(count)=y;
end
PP=Pfun(YY);
f=target(TT,PP);
if isplot
figure;
plot(TT,PP,'r');
xlim([min(TT),max(TT)]);
ylim([min(PP),max(PP)]);
end
end
19.V_in_2.m
function V=V_in_2(w,t)
A=2.5^2*pi;
V0=20+2.413*A*2;
rho=-2.413*cos(w*t)+2.413;
V=V0-rho*A;
end
20.yfun_1.m
function [y,m]=yfun_1(t,y,m,w,dt)
Q_in=@Q_in_1;
Q_out=@Q_out_1;
V_0=500*(pi/4)*(10)^2;
y=y+((1/V_0)*(0.87*Q_in(w,t,Pfun(y))-y*Q_out(t,Pfun(y))))*dt;
end
21.yfun_2.m
function [y,m]=yfun_2(t,y,m,w,dt)
if (mod(t,2*pi/w)<dt)
m=V_in_2(w,0)*0.850;
end
V_0=500*(pi/4)*(10)^2;
V_A=V_in_2(w,t);
rho_A=m/V_A;
m_out=Q_out_2(t,y)*y;
P_A=Pfun(rho_A);
P_N=Pfun(y);
if (P_A>P_N)
Q_A=0.85*0.7*0.7*3.1415*sqrt(2*(P_A-P_N)/rho_A);
y=y+(1/V_0)*(rho_A*Q_A-m_out)*dt;
m=m-rho_A*Q_A*dt;
if (m<0)
m=0;
end
else
if (rho_A==0)
pause;
end
y=y+(1/V_0)*(-m_out)*dt;
end
end
22.yfun_3.m
function [y,m]=yfun_3(t,y,m,w,dt)
if (mod(t,2*pi/w)<dt)
m=V_in_2(w,0)*0.850;
end
V_0=500*(pi/4)*(10)^2;
V_A=V_in_2(w,t);
rho_A=m/V_A;
m_out=Q_out_2(t,y)*y*2;
P_A=Pfun(rho_A);
P_N=Pfun(y);
if (P_A>P_N)
Q_A=0.85*0.7*0.7*3.1415*sqrt(2*(P_A-P_N)/rho_A);
y=y+(1/V_0)*(rho_A*Q_A-m_out)*dt;
m=m-rho_A*Q_A*dt;
if (m<0)
m=0;
end
else
if (rho_A==0)
pause;
end
y=y+(1/V_0)*(-m_out)*dt;
end
end
23.yfun_4.m
function [y,m]=yfun_4(t,y,m,P_L,dt)
global flag;
w=0.05045;
if (mod(t,2*pi/w)<dt)
m=V_in_2(w,0)*0.850;
end
V_0=500*(pi/4)*(10)^2;%高压油管容积
V_A=V_in_2(w,t);%高压油泵内体积
rho_A=m/V_A;%高压油泵内密度
m_out_1=0;%喷油嘴堵塞
m_out_2=Q_out_3(Pfun(y),y)*y;%减压阀出油质量
m_out_3=Q_out_2(t,y)*y*2;%喷油嘴修复后出油质量
P_A=Pfun(rho_A);%高压油泵内压力
P_N=Pfun(y);%高压油管内压力
P_H=180;
if(P_N>P_H)
flag=1;
end
if (flag==0)%喷油嘴堵塞
if (P_A>P_N)
Q_A=0.85*0.7*0.7*pi*sqrt(2*(P_A-P_N)/rho_A);
P_A=Pfun(rho_A);
y=y+(1/V_0)*(rho_A*Q_A-m_out_1)*dt;
P_N=Pfun(y);
m=m-rho_A*Q_A*dt;
if (m<0)
m=0;
end
else
y=y+(1/V_0)*(-m_out_1)*dt;
end
end
%喷油嘴和减压阀同时打开
if (flag==1&&P_N>=P_L)
if (P_A>P_N)
Q_A=0.85*0.7*0.7*pi*sqrt(2*(P_A-P_N)/rho_A);
y=y+(1/V_0)*(rho_A*Q_A-m_out_2-m_out_3)*dt;
m=m-rho_A*Q_A*dt;
if (m<0)
m=0;
end
else
y=y+(1/V_0)*(-m_out_2-m_out_3)*dt;
end
end
%减压阀关闭
if (P_N<P_L)
if (P_A>P_N)
Q_A=0.85*0.7*0.7*pi*sqrt(2*(P_A-P_N)/rho_A);
y=y+(1/V_0)*(rho_A*Q_A-m_out_3)*dt;
m=m-rho_A*Q_A*dt;
if (m<0)
m=0;
end
else
y=y+(1/V_0)*(-m_out_3)*dt;
end
end
end