国赛2019A——高压油管问题

问题背景

 

 

各参数计算公式

弹性模量计算压力

function k = E2P(P,E)
    k = polyfit(P,E,2);
end

进油公式

function Q = Q_in(P1,P2,rho,t,T)
    %% 仿真时间t时刻的进油量流量函数
    C = 0.85;
    A = pi*1.4*1.4/4;
    t_convert = mod(t,T+10);
    if t_convert <= T
        Q = C*A*sqrt(2*(P1-P2)/rho);
    else
        Q = 0;
    end
end

出油公式

function Q = Q_out(t)
    %% 仿真时间t时刻的出油量流量函数
    % 假设喷油嘴每秒以均匀频率工作10次，此时喷油量可看作周期函数
    t_convert = mod(t,100); % 周期转换
    if t_convert <= 0.2 && t_convert >= 0
        Q = 100*t_convert;
    elseif 0.2 <t_convert && t_convert < 2.2
        Q = 20;
    elseif 2.2 <= t_convert && t_convert <= 2.4
        Q = 240-100*t_convert;
    else
        Q = 0;
    end
end

欧拉法计算密度

function rho = rho_iter(rho0,dP,E,P)
    if P > 100
        rho = rho0*(1+dP/E);
    else
        rho = rho0*(1-dP/E);
    end
end

密度计算压力

clear;clc
data = readmatrix('附件3-弹性模量与压力.xlsx','Range','A2:B402');
P = data(:,1);
E = data(:,2);
rho = zeros(size(P));
ind = find(P==100);
rho(ind) = 0.85;
n = length(P);
dP = mean(diff(P,1));
for i = ind:n-1
    rho(i+1) = rho_iter(rho(i),dP,E(i),P(i));
end
for i = ind:-1:2
    rho(i-1) = rho_iter(rho(i),dP,E(i),P(i));
end
k = polyfit(rho,P,2);
P_test = k(1)*rho.^2+k(2)*rho+k(3);
plot(rho,P)
hold on 
plot(rho,P_test)

function P = rho2P(rho)
    k = [10952.7722859703,-15955.0071828209,5749.08734491654];
    P = k(1)*rho.^2+k(2)*rho+k(3);
end

压力计算

function [P,rho,M] = P_vessel(P1,P2,t,T,rho1,rho2,M,V,Qout,dt)
    Qin = Q_in(P1,P2,rho1,t,T);
    M = M+dt*rho1*Qin-dt*rho2*Qout;
    rho = M/V;
    P = rho2P(rho);
end

各时刻压力仿真计算

function [t,P] = pipe_sim(P0,T,t_sim)
    %% 基本参数设置
    % 油管参数
    L = 500; % 油管长度
    d_pipe = 10; % 油管直径
    V = L*pi*(d_pipe/2)^2; % 油管体积
    % 仿真时间参数
    dt = 0.1;
    t = 0:dt:t_sim; t = t';
    % 状态量参数
    rho0 = 0.8707;
    Qout = zeros(size(t)); % 出油量
    P = zeros(size(t)); P(1) = 100; % 管内压力
    rho = zeros(size(t)); rho(1) = 0.85; % 管内燃油密度
    M = zeros(size(t)); M(1) = rho(1)*V; % 管内燃油质量
    for i = 1:length(t)
        Qout(i) = Q_out(i);
    end
    for i = 1:length(t)-1
        [P(i+1),rho(i+1),M(i+1)] = P_vessel(P0,P(i),t(i),T,rho0,rho(i),M(i),V,Qout(i),dt);
    end 
end

粒子群法计算单次阀门开启时间

function [e,T] = scope_T(P_target,t_sim,P0)
    %% 粒子群算法中的预设参数
    n = 10; % 粒子数量
    narvs = 1; % 变量个数
    c1 = 2;  % 每个粒子的个体学习因子，也称为个体加速常数
    c2 = 2;  % 每个粒子的社会学习因子，也称为社会加速常数
    w = 0.9;  % 惯性权重
    K = 50;  % 迭代的次数
    vmax = 1.2; % 粒子的最大速度
    x_lb = 0; % x的下界
    x_ub = 1; % x的上界

    %% 初始化粒子的位置和速度
    x = zeros(n,narvs);
    for m = 1: narvs
        x(:,m) = x_lb(m) + (x_ub(m)-x_lb(m))*rand(n,1);   
    end
    v = -vmax + 2*vmax .* rand(n,narvs);  

    %% 计算适应度
    fit = zeros(n,1); 
    for m = 1:n  
        fit(m) = Obj_fun(P0,x(m,:),t_sim,P_target);  
    end
    pbest = x;  
    ind = find(fit == max(fit), 1); 
    gbest = x(ind,:);  

    %% 迭代K次来更新速度与位置
    fitnessbest = ones(K,1);
    for d = 1:K  
        for m = 1:n 
            v(m,:) = w*v(m,:) + c1*rand(1)*(pbest(m,:) - x(m,:)) + c2*rand(1)*(gbest - x(m,:));  % 更新第i个粒子的速度
            % 如果粒子的速度超过了最大速度限制，就对其进行调整
            for j = 1: narvs
                if v(m,j) < -vmax(j)
                    v(m,j) = -vmax(j);
                elseif v(m,j) > vmax(j)
                    v(m,j) = vmax(j);
                end
            end
            x(m,:) = x(m,:) + v(m,:); % 更新第i个粒子的位置
            % 如果粒子的位置超出了定义域，就对其进行调整
            for j = 1: narvs
                if x(m,j) < x_lb(j)
                    x(m,j) = x_lb(j);
                elseif x(m,j) > x_ub(j)
                    x(m,j) = x_ub(j);
                end
            end
            fit(m) = Obj_fun(P0,x(m,:),t_sim,P_target);
            if fit(m) > Obj_fun(P0,pbest(m,:),t_sim,P_target) 
                pbest(m,:) = x(m,:);  
            end
            if  fit(m) > Obj_fun(P0,gbest,t_sim,P_target) 
                gbest = pbest(m,:); 
            end
        end
        fitnessbest(d) = Obj_fun(P0,gbest,t_sim,P_target); 
    end
    T = gbest;
    e = fitnessbest;   
end


function e = Obj_fun(P0,x,t_sim,P_target)
    [~,P] = pipe_sim(P0,x,t_sim);
    P_mean = mean(P(80000:100000));
    e = -abs(P_mean-P_target);
end

仿真过程

clear;clc

%% 仿真过程
t_sim = 10000;
P0 = 160;
P_target = input('请输入稳定时的压强：');
[e,T] = scope_T(P_target,t_sim,P0); e = -e;
[t,P] = pipe_sim(P0,T,t_sim);

%% 作图
figure(1)
plot(t,P,'k')
xlabel('时间/ms')
ylabel('压力/MPa')
axis([0 t_sim 0 200])
figure(2)
plot(e,'k')
xlabel('迭代次数')
ylabel('误差')