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

原创已于 2025-12-31 14:37:15 修改 · 513 阅读

10 ·

本内容遵循CC 4.0 BY-SA版权协议

GEO检测

标签

#matlab #粒子群 #微分方程 #二项式拟合 #高压油管

于 2023-07-02 20:30:22 首次发布

该文详细阐述了一种用于计算管道内部压力的模拟方法。通过弹性模量计算压力、进油和出油的流量函数，以及使用欧拉法计算密度，建立了压力变化模型。此外，利用粒子群优化算法寻找阀门最佳开启时间，以达到目标压力。最终，通过仿真展示了系统在不同时间的压力变化和算法的优化效果。

问题背景

高压油管的压力控制是一道自控方向的题目。通过调节燃油压力和流量，确保喷油器精准喷油，同时避免压力波动对系统造成损害。通过对不同喷油策略的压力时序仿真，得到最优的单向阀开关时间。

各参数计算公式

弹性模量计算压力

通过二项式拟合，求解弹性模量和压力间的非线性关系。

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

欧拉法计算密度

欧拉法是一种数值解常微分方程的经典方法，由数学家欧拉提出。其核心思想是用离散的差分近似代替连续的微分，通过迭代逐步逼近真实解。对于初值问题 $\frac{dy}{dt} = f(t, y)$ ，欧拉法的迭代公式为： $y_{n+1} = y_n + h \cdot f(t_n, y_n)$ ，其中 ( h ) 为步长。根据数值分析中的欧拉法，计算密度的迭代变化。得到密度迭代函数。

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

各时刻压力仿真计算

按时间维度对压力值进行仿真。for循环压力值函数。得到对应的稳压值。

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

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

粒子群算法（Particle Swarm Optimization, PSO）是一种基于群体智能的优化算法，模拟鸟群或鱼群的社会行为。算法通过个体间的协作与信息共享寻找最优解，适用于连续空间优化问题。

每个粒子代表解空间中的一个候选解，具有位置和速度两个属性。粒子根据自身历史最优解和群体历史最优解动态调整飞行方向与速度，逐步逼近全局最优解。其核心公式为速度更新公式：

$v_{i}(t+1) = w \cdot v_{i}(t) + c_1 \cdot r_1 \cdot (pbest_i - x_i(t)) + c_2 \cdot r_2 \cdot (gbest - x_i(t))$

其中 $w$ 为惯性权重， $c_1c_2$ 为学习因子， $r_1r_2$ 为随机数。

PSO算法参数少、收敛快，广泛应用于函数优化、神经网络训练、控制系统设计等领域。其变种包括带惯性权重的PSO、混合PSO等，以平衡全局探索与局部开发能力。

在给定初压值、稳压值与仿真时间，通过粒子群算法对最优（min）的阀门开启时间T进行寻优。寻优目标为稳压值与给定稳压值间距离的绝对值。（系泊系统那篇文章也用过）

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('误差')

得到的仿真图如下：

第一张和第二张分别是T=5s和10s达到150MPa稳压的情况，第三张和第四张为通过寻优确定的达到100MPa和150MPa稳压所需的最短阀门开启时间T。