Finite Element Model Updating for Bridges Based on Adaptive Nested Sampling and Bayesian Theory
-
摘要:
在基于有限元模型的桥梁健康监测中,贝叶斯模型修正技术通常被用于量化有限元模型中重要参数的不确定性,以解决模型修正中由于测量误差、建模误差、计算误差等造成的非唯一解问题. 为解决由于大量调用有限元模拟运算,导致修正效率低下的问题,基于自适应嵌套抽样(ANS)算法,提出一种贝叶斯模型修正方法. 该方法利用模态参数构建概率目标函数,并采用 ANS 算法对其进行逼近,ANS 保留了嵌套抽样(NS)的性质,通过逐层缩小抽样范围,使得样本最终逼近最优参数;通过逐层近似,将高维积分问题转化为简单的一维积分问题,简化了证据值和后验概率密度值的计算过程;在此基础上,ANS 算法在迭代过程中通过自适应地调整样本数量,减少对有限元模型的调用;最后,对一座人行桁架桥进行了贝叶斯有限元模型修正试验. 结果表明:在相同算法参数设置下,ANS 算法相比传统 NS 算法降低了约 84% 的有限元模拟调用次数,节省了约 86% 计算时间,并能获得同等精度的不确定性修正结果.
Abstract:In bridge health monitoring based on finite element models, Bayesian model updating techniques are commonly used to quantify the uncertainties of important parameters in the finite element models, so as to address the issue of non-uniqueness in model updating caused by measurement errors, modeling errors, computational errors, etc. To resolve the problem of low efficiency in model updating due to the large number of finite element simulations required, a Bayesian model updating method based on an adaptive nested sampling (ANS) algorithm was proposed. The method used the modal parameters to construct the probability objective function and adopted the ANS algorithm to approximate it. ANS retained the nature of nested sampling (NS), which made the samples ultimately approximate the optimal parameters by narrowing the sampling range layer by layer, and it simplified the computation process of the evidence value and the a posteriori probability density value by transforming the high-dimensional integration problem into a simple one-dimensional integration problem through layer-by-layer approximation. On this basis, the ANS algorithm could also reduce the call of the finite element model by adaptively adjusting the number of samples during the iteration process. Finally, a pedestrian truss bridge was used as a case study for Bayesian finite element model updating experiments. The results demonstrate that under the same algorithm parameter settings, the ANS algorithm reduces the number of finite element simulation calls by approximately 84% compared to the traditional NS algorithm. This leads to approximately 86% computational time savings while obtaining uncertainty updating results with equal accuracy.
-
空车调运是实现铁路资源优化配置的重要内容,空车调运的优劣直接影响铁路车辆的利用效率、市场需求的满足程度以及运输组织的现代化水平.
近年来空车调运优化问题的研究多从供需及环境的随机性出发建立动态随机规划模型. 文献[1]基于广义旅速建立时空服务网络,提出了空车调配多阶段动态优化方法;在此基础上,文献[2]着重考虑了实际运输生产中的能力约束,建立了考虑车种替代动态规划模型;文献[3]则将空车调配分为基于固定需求的优化调整和基于客户实时需求的策略再优化两个阶段,解决了空车需求的动态变化问题;文献[4]基于供需的不确定性,为约束条件和目标函数设置了一定的置信水平,建立了多车种调配优化模型;文献[5]考虑了路网中转能力和通行能力不确定性的影响,建立多目标随机期望约束模型;文献[6]则针对列车走行时间的不确定性,提出了鲁棒优化模型,并利用对称和对偶变化降低了模型求解难度,提高了求解效率.
在上述研究的基础上,考虑到实际铁路生产中空车的产生以及空车的需求是一个时空问题,除了物理层面的连接外,还应当结合空车与列车接续时间关系进行列车之间的空车配流. 同时,网络中列车的旅行时间以及车站的技术作业时间具有不确定性,时间的波动会直接对供需网络的接续关系造成影响. 因此,本文基于鲁棒优化理论对多车种、可替代条件下的空车调运问题进行了研究.
1. 问题描述与建模
1.1 问题描述
供应站与需求站的车站技术作业时间以及站间旅行时间不确定性时,如何确定需求站的空车来源、类型与空车数量是鲁棒接续时间要求下空车调运的核心问题. 模型中包含下述参数.
定义集合及索引:I、i分别为空车供应站集合和索引,
i∈I ;J、j分别为空车需求站集合及索引,j∈J ;E、e分别为供应站内产生空车的到达列车集合及索引,e∈E ;F、f分别为供应站发出的可挂运空车的列车集合及索引,f∈F ;G 、g 分别为需求站内需使用到达空车装车的发出列车集合及索引,g∈G ;η 为实际的空车类型,χ 为考虑车种替代后按用途归类的空车类型,空车类型集合H包括平车(N)、棚车(P)和敞车(C),H={N,P,C} ,且η,χ∈H .定义供应站相关变量:
ci 、ti 分别为供应站i的等待车小时费用及车站技术作业时间;Sie,η 、Sif 分别为供应站i内第e列到达列车产生η 型空车的数目及第f 列发出列车的最大空车挂运数目;Aie 、Lif 分别为供应站i第e列到达列车的到达时刻及第f列发出列车的最晚编组时刻.定义需求站相关变量:pj、cj和tj分别为需求站j的空车装车后产生收益、空车等待车小时费用及车站技术作业时间;Sjg,χ为空车需求站j第g列发出列车需要的
χ 型空车数量;Ljg为空车需求站j第g列发出列车的最晚编组时刻.定义0-1变量: xi,e,f表示供应站i的第e列到达列车与供应站i的第f列发出列车间是否存在接续关系,xi,e,f = 1为是,xi,e,f = 0为否;yif,jg = 0表示供应站i的第f列发出列车与需求站j的第g列发出列车间是否存在接续关系,yif,jg = 1为是,yif,jg = 0为否.
定义决策变量:di,e,f,η为供应站i的第e列到达列车为供应站i发出的第f列列车供应的η型空车的数量;dif,jg,χ为供应站i第f列发出列车为需求站j发出的第g列列车供应的χ型空车的数量;dif,jg,η为供应站i第f列发出列车为需求站j第g列发出列车供应的η型空车的数量.
另外定义了空车供应站i和空车需求站j间的空车输送费用cij及站间旅行时间tij.
1.2 名义模型构建
名义模型是将供应站与需求站车站技术作业时间以及站间旅行时间看作确定值时的空车调运模型.
1) 目标函数
以空车调运效益最大化为目标函数,如式(1)所示,其中:等号右端第1项表示空车到达需求站后由需求站装车发出所产生收益;第2项表示空车调运时产生的路径消耗费用[7];第3项表示空车在需求站的等待费用;第4项表示空车在供应站的等待费用.
C1=max{∑χ∑i∑j∑f∑gyif,jgdif,jg,χpj−∑χ∑i∑j∑f∑gyif,jgdif,jg,χcij−∑χ∑i∑j∑f∑gyif,jgdif,jg,χcj[Ljg−(Lif+tij+tj)]−∑η∑i∑e∑fxi,e,fdi,e,f,ηci(Lif−ti−Aie)}. (1) 2) 约束条件
式(2)表示供应站i的第e列到达列车给该站各发出列车f提供的η型空车总数不得超过第e列到达列车产生η型空车总数.
∑fxi,e,fdi,e,f,η⩽Sie,η. (2) 式(3)表示供应站i的各次到达列车e为该站第f列发出列车供应的空车总数不得超过第f列发出列车的最大空车挂运数目.
∑η∑exi,e,fdi,e,f,η⩽Sif. (3) 式(4)表示
dif,jg,η 中η 型箱的数量等于其实际替换的各车型dif,jg,χ 的数量总和. 设Rη 为η 型箱可替代的车型集合,参考文献[8],平车代替敞车有90%概率,棚车代替平车和敞车有75%的可能,敞车代替棚车有30%的可能,故有RN = {N,C},RP = {N,P,C},RC = {P,C}.dif,jg,η=∑χ∈Rηdif,jg,χ. (4) 式(5)表示供应站i发出的第f列列车中挂运的各型空车均用于满足需求站的空车需求[9-10].
∑exi,e,fdi,e,f,η=∑gyif,jgdif,jg,η. (5) 式(6)表示各供应站发出的列车为需求站j发出的第g列列车提供的χ型空车数量等于该列车对χ型空车的需求量.
∑i∑fyif,jgdif,jg,χ = Sjg,χ. (6) 式(7)表示列车间空车输送数量非负.
{di,e,f,η⩾0,dif,jg,η⩾0. (7) 式(8)表示如果供应站i的第e列到达列车在供应站i完成相关技术作业后,到达调车场的时刻不晚于i站第f列发出列车的最晚编组时刻,则该组供应站到达列车与发出列车之间满足最小接续时间要求,即xi,e,f = 1,否则,xi,e,f = 0. 供应站到达列车产生的空车来自重车卸车时,ti = (ti,dj + ti,zx) + (ti,jt + ti,zx) + (ti,xc + ti,zx);到达列车本身挂有空车时,ti = (ti,dj + ti,zx) + (ti,jt + ti,zx);其中:ti,dj 、 ti,zx 、ti,jt和ti,xc分别为i站的到达技术作业时间标准、转线作业时间标准、解体作业时间标准(含推峰时间)及卸车作业时间标准.
xi,e,f={1,Aie+ti⩽Lif,0,Aie+ti>Lif. (8) 式(9)表示若供应站i发出的第f列列车在需求站j完成相关技术作业后,到达调车场的时刻不晚于j站发出第g列列车的最晚编组时刻,即yie,jg = 1,否则,yie,jg = 0. 列车在需求站j的技术作业时间tj = (tj,dj + tj,zx) + (tj,jt + tj,zx) + (tj,zc + tj,zx),其中:tj,dj、tj,zx、tj,jt、tj,zc分别为j站到达技术作业时间标准、转线作业时间标准、解体作业时间标准及装车作业的时间标准.
yif,jg={1,Lif+tij+tj⩽Ljg,0,Lif+tij+tj>Ljg. (9) 式(10)表示只有当供应站i第f列发出列车存在满足接续时间要求的空车来源时,即xi,e,f = 1时,才考虑该发出列车与需求站发出列车的接续关系.
yif,jg⩽∑exi,e,f. (10) 为便于说明,本文定义了属于关系、连接关系和车流关系,如图1所示,图中:三角框为供应站,正方形框为需求站,圆圈为列车,后同.
定义1 属于关系指供应站到达列车e、发出列车f与供应站i之间的关系以及需求站发出列车g与需求站j之间的关系.
定义2 连接关系为供应站发出的列车f与该列车到达的需求站j之间的关系,供应站发出的任一可挂运空车的列车有且只与一个需求站建立连接关系.
定义3 车流关系是指当某一列车中的车辆可以成为另一列车中的车流来源时,两列列车之间存在的关系. 车流关系可以定义为两部分:
1) 供应站到达列车e与该站发出列车f之间的车流关系定义为第1阶段车流关系;
2) 供应站发出的列车f与需求站发出的列车e之间的车流关系定义为第2阶段车流关系.
车流关系由约束条件式(8)~(10)确定.
1.3 鲁棒优化模型构建
实际运输中由于设备因素、环境因素和人为因素的影响,列车站间旅行时间以及车站技术作业时间是不确定的,若以固定的旅行时间和技术作业时间作为空车接续时间的判断依据,当接续条件
Lif − (Aie+ti)⩾0 或Ljg − (Lif+tij+tj)⩾0 ,但取值较小时,技术作业时间和旅行时间的波动可能使得既有的接续关系失效,导致求得方案不可行. 因此,基于Dimitris等[11]提出的鲁棒优化理论建立了空车调运鲁棒优化模型.为描述站间旅行时间和车站技术作业时间的不确定性,令
˜tij∈[tij,tij+ˆti] 、˜ti∈[ti,ti+ˆti] 、˜tj∈[tj,tj+ˆtj] ,ˆtij 、ˆti 、ˆtj 分别为tij、ti、tj的时间最大波动量.引入
αij=(˜tij−tij)/ˆtij 、βi=(˜ti−ti)/ˆti 和γj=(˜tj−tj)/ˆtj 分别表示tij 、ti 和tj 波动程度的系数,αij,βi,γj∈[0,1] . 令参数集合为U、V、W,αij ∈U ,βi ∈V ,γj ∈W ,可得{˜tij=tij+αijˆtij,˜ti=ti+βiˆti,˜tj=tj+γjˆtj. (11) 若
˜tij、˜ti、˜tj 均波动,则模型过于保守,为避免此情形,在模型中引入Γ1 、Γ2 和Γ3 来调整模型的鲁棒性,以此作为波动下限值,参数的不确定集为∑i∑jαij⩾Γ1 ,∑iβi⩾Γ2 ,∑jγj⩾Γ3 .Γ1 、Γ2 和Γ3 保证了各项时间均存在一定的不确定性,且取值越大,时间的不确定程度越大,调整Γ1 、Γ2 和Γ3 的取值可以避免模型过度鲁棒. 构建的鲁棒优化模型如下.1) 目标函数
鲁棒优化模型的目标函数如式(12),约束关系如式(13)所示,结合波动程度与波动下限约束,式(12)中的C2如(14)所示,矩阵表达如式(15)所示.
max{C1+C2}, (12) s.t.{αij + δij=1,βi+δi=1,γj+δj=1,∑i∑j(θij−αij)=−Γ1,∑i(θi−βi)=−Γ2,∑j(θj−γj)=−Γ3, (13) C2=min{∑χ∑i∑j∑f∑gyif,jgdif,jg,χci(αijˆtij+γjˆtj)+∑η∑i∑j∑fxi,e,fdi,e,f,ηβiˆti}, (14) \left[\begin{array}{ccccccccc}{{\boldsymbol{E}}}_{nm}& 0& 0& {{\boldsymbol{E}}}_{nm}& 0& 0& 0& 0& 0\\ 0& {{\boldsymbol{E}}}_{n}& 0& 0& {{\boldsymbol{E}}}_{n}& 0& 0& 0& 0\\ 0& 0& {{\boldsymbol{E}}}_{m}& 0& 0& {{\boldsymbol{E}}}_{m}& 0& 0& 0\\ -{{\boldsymbol{Z}}}_{nm}& 0& 0& 0& 0& 0& {{\boldsymbol{Z}}}_{nm}& 0& 0\\ 0& -{{\boldsymbol{Z}}}_{n}& 0& 0& 0& 0& 0& {{\boldsymbol{Z}}}_{n}& 0\\ 0& 0& -{{\boldsymbol{Z}}}_{m}& 0& 0& 0& 0& 0& {{\boldsymbol{Z}}}_{m}\end{array}\right]\left[\begin{array}{c}{{\boldsymbol{U}}}_{nm}\\ {{\boldsymbol{V}}}_{n}\\ {{\boldsymbol{W}}}_{m}\\ {{\boldsymbol{\delta}} }_{nm}\\ {{\boldsymbol{\delta}} }_{n}\\ {{\boldsymbol{\delta}} }_{m}\\ {{\boldsymbol{\theta}} }_{nm}\\ {{\boldsymbol{\theta}} }_{n}\\ {{\boldsymbol{\theta}} }_{m}\end{array}\right]=\left[\begin{array}{l}{{\boldsymbol{Z}}}_{nm}^{\text{T}}\\ {{\boldsymbol{Z}}}_{n}^{\text{T}}\\ {{\boldsymbol{Z}}}_{m}^{\text{T}}\\ -{\varGamma }_{1}\\ -{\varGamma }_{2}\\ -{\varGamma }_{3}\end{array}\right], (15) 式(13)、(15)中:
{\delta _{ij}} 、{\delta _i} 、{\delta _j} 、{\theta _{ij}} 、{\theta _i} 、{\theta _j} 为松弛系数;n 为供应站数目;m 为需求站数目;{{\boldsymbol{E}}_n} 为n \times n 的单位矩阵;{{\boldsymbol{Z}}_n} 为元素全为1的1 \times n 矩阵;{{\boldsymbol{\delta}} _n} = {({\delta _1},{\delta _2}, \cdots ,{\delta _n})^{\rm T}} ;{{\boldsymbol{\theta}} _n} = {({\theta _1},{\theta _2}, \cdots ,{\theta _n})^{\rm T}} ;{{\boldsymbol{U}}_{nm}} = (\alpha _{11},\alpha_{12}, \cdots , \alpha _{nm})^{\rm T} ;{{\boldsymbol{V}}_n} = {(\beta_1,\beta _2, \cdots ,{\beta _n})^{\rm T}} ;{{\boldsymbol{W}}_m} = {(\gamma_1,{\gamma _2}, \cdots ,{\gamma _m})^{\rm T}} .证明一个矩阵为完全幺模矩阵的充分条件为:① 矩阵中元素仅包含0、1、−1;② 矩阵每列最多两个非0元素;③ 矩阵的行可分划成两个子集,使得同列中两个非零元素符号相同时,对应的两行在不同的行子集中,当符号不同时,对应的两行在同一行子集中. 易得:约束矩阵式(15)满足条件 ①、②,并且当矩阵行划分为(1,4)和(2,3,5,6)两个行子集时满足条件 ③,故该约束矩阵为完全幺模矩阵. 结合完全幺模矩阵的性质,当
{\varGamma _1} 、{\varGamma _2} 和{\varGamma _3} 均取整数时,式(13)、(14)的最优解{\alpha ^{\text{*}}} 、{\beta ^{\text{*}}} 、{\lambda ^{\text{*}}} 也为整数,因{\alpha _{ij}},\;{\beta _i},\;{\gamma _j}\geqslant 0 ,且{\alpha _{ij}},\;{\beta _i}, {\gamma _j} \leqslant 1 ,故其最优值取0或1.2) 约束条件
原约束条件中式(2)~(7)、(10)不变,式(8)、(9)分别变更为式(16)、(17).
{x_{i,e,f}} = \left\{ \begin{gathered} 1,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {A_{ie}} + {t_i} + {\beta _i}\hat{t} {{_i}} \leqslant {L_{if}}, \hfill \\ 0,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {A_{ie}} + {t_i} + {\beta _i}\hat{t} {{_i}} > {L_{if}}, \hfill \\ \end{gathered} \right. (16) {y_{if,jg}} = \left\{ \begin{gathered} 1,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {L_{if}} + {t_{ij}} + {a_{ij}}\hat{t} {{_{ij}}} + {t_j} + {\gamma _j}\hat{t} {{_j}} \leqslant {L_{jg}}, \hfill \\ 0,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {L_{if}} + {t_{ij}} + {a_{ij}}\hat{t} {{_{ij}}} + {t_j} + {\gamma _j}\hat{t} {{_j}} > {L_{jg}}{\kern 1pt} . \hfill \\ \end{gathered} \right. (17) 2. 求解算法
2.1 名义模型求解算法
名义模型为非线性整数规划模型,求解困难,因此,求解前结合约束关系式(8) ~ (10)对模型进行简化.
以算例中的供应站3为例,为便于表示,此处仅保留三列供应站发出列车,变化前后网络结构如图2所示. 与原网络相比,新网络中仅保留了
{x_{i,e,f}} = 1 和{y_{ie,jg}} = 1 时所对应的决策变量{d_{i,e,f,\eta }} 和{d_{if,jg,\eta }} ,进而去除模型中的约束条件式(8)~(10),转化为线性规划模型,可调用CPLEX进行求解.2.2 鲁棒优化模型求解算法
鲁棒优化模型中引入了变量
{\alpha _{ij}} 、{\beta _i} 、{\gamma _j} ,这些变量会导致供需网络中既有车流关系的变化,但对于每组{\alpha _{ij}} 、{\beta _i} 、{\gamma _j} 存在唯一的网络结构,求解时可以看作名义模型.遍历满足约束的每组
{\alpha _{ij}} 、{\beta _i} 、{\gamma _j} 是最准确的求解方法,但当{\varGamma _1} 、{\varGamma _2} 、{\varGamma _3} 取值较高时模型的求解复杂度会剧增,导致求解困难. 研究发现,只有当网络结构中的车流关系减少时,最优调运方案才可能改变,故设计下述步骤求解. 求解时{\varGamma _1} 、{\varGamma _2} 、{\varGamma _3} 和\hat{t} {{_{ij}}} 、\hat{t} {{_i}} 、\hat{t} {{_j}} 为已知值.步骤1 结合式(8)~(10),得出固定出行时间下供需网络中存在的车流关系,令供应站i的第1阶段车流关系集合为
{{\varOmega} _{{\kern 1pt} {\kern 1pt} 1,i}} ,供应站i的第2阶段车流关系集合为{{\varOmega} _{{\kern 1pt} {\kern 1pt} {\text{2}},i}} ,形式如下:{{\varOmega} _{{\kern 1pt} {\kern 1pt} 1,i}} = [(1,1),(1,2), \cdots ,(e,f)] , (18) \begin{split} {{\varOmega} _{{\kern 1pt} {\kern 1pt} {\text{2,}}i}} = [(1,1,1),(1,1,2), \cdots ,(f,j,g)], \end{split} (19) 式中:Ω1,i中(
e,f )表示i站第e列到达列车可以为i站第f列发出列车供应空车;Ω2,i中(f,j,g)表示i站发出的第f列列车可以为j站发出的第g列列车供应空车.步骤2 令集合
U 、V 、W 中的各项元素均为1,即所有的{t_{ij}} 、{t_i} 、{t_j} 都存在波动,结合式(10)、(16)、(17),得出绝对鲁棒下的网络中存在的车流关系{{\varPi} _{{\kern 1pt} {\kern 1pt} 1,i}} 、{{\varPi} _{{\kern 1pt} {\kern 1pt} {\text{2}},i}} ,形式同{{\varOmega} _{{\kern 1pt} {\kern 1pt} 1,i}} 、{{\varOmega} _{{\kern 1pt} {\kern 1pt} {\text{2}},i}} .步骤3 比较
{{\varPi} _{{\kern 1pt} {\kern 1pt} 1,i}} 、{{\varPi} _{{\kern 1pt} {\kern 1pt} {\text{2}},i}} 和{{\varOmega} _{{\kern 1pt} {\kern 1pt} 1,i}} 、{{\varOmega} _{{\kern 1pt} {\kern 1pt} {\text{2}},i}} 中的车流关系数量,确定此时减少的车流关系. 若{{\varPi} _{{\kern 1pt} {\kern 1pt} 1,i}} \subsetneqq {{\varOmega} _{{\text{1,}}i}} ,说明第1阶段车流关系发生了变化,若{{\varPi} _{{\kern 1pt} {\kern 1pt} {\text{2}},i}} \subsetneqq {{\varOmega} _{{\kern 1pt} {\kern 1pt} {\text{2}},i}} ,说明第2阶段车流关系发生变化. 假设此时减少的车流关系为{{\varOmega} _{{\kern 1pt} {\kern 1pt} 1,1}} 的(1,1),{{\varOmega} _{{\text{1,2}}}} 的(1,2),{{\varOmega} _{{\kern 1pt} {\kern 1pt} {\text{2}},2}} 的(1,1,1),{{\varOmega} _{{\kern 1pt} {\kern 1pt} {\text{2}},{\text{3}}}} 的(2,2,1).步骤4 找出减少的车流关系中包含的供应站集合I1、需求站集合J1以及连接关系集合O,以步骤3中假设为例,减少的车流关系中包含的供应站为I1 = {1,2,3},包含的需求站为J1 = {1,2},包含的连接关系为
{O} = \{ (2,1), (3,2)\} .步骤5 令集合
U 、V 、W 中的各项元素均为0. 同时,根据供应站{n_1} 、需求站{n_2} 、连接关系{n_3} 及波动下限{{\varGamma _1}、{\varGamma _2}、{\varGamma _3}} ,按下述方法确定{\beta _i} 、{\gamma _j} 、{\alpha _{ij}} 的取值:1) 若
{\varGamma _1} \geqslant \left| {{n_1}} \right| (\left| {{n_1}} \right| 表示{n_1} 中元素个数,其余同),则令{I_1} 内各供应站的系数{\beta _i} 为1,以{I_1}{\rm{ = }}\left\{ {{\rm{1,2,3}}} \right\} 为例即{\beta _1} = 1 、{\beta _2} = 1 、{\beta _3} = 1 ;若{\varGamma _1} < \left| {{I_1}} \right| 则得到从{I_1} 中取{\varGamma _1} 个元素的全排列集合{N_1} ,令每组排列的系数为1;2) 针对
{\varGamma _2} 的调整方法同T1,当{\varGamma _{\text{2}}} < \left| {{J_{\text{1}}}} \right| 时,得到从{J_{\text{1}}} 中取{\varGamma _{\text{2}}} 个元素的全排列集合{N_{\text{2}}} ;3) 若
{\varGamma _{\text{3}}} \geqslant \left| {{O}} \right| ,则令{O} 中各元素的系数{\alpha _{ij}} = 1 ,以步骤3中假设为例,{\alpha _{21}} = 1 ,{\alpha _{32}} = 1 ;否则,得到从{O} 中取{\varGamma _3} 个元素的全排列组合{N_3} ,令每组排列的系数为1.步骤6 输入由
{N_1} 、{N_2} 、{N_3} 中元素组成的{\alpha _{ij}} 、{\beta _i} 、{\gamma _j} 排列,调用CPLEX求解各网络结构下的最佳效益值,取其中的最小值作为鲁棒解.该算法在给定波动量
\hat{t} {{_{ij}}} 、\hat{t} {{_i}} 、\hat{t} {{_j}} 的情况下,找出最劣网络结构较理想网络结构中减少的车流关系,以此为入手点找出波动下限{\varGamma _1} 、{\varGamma _2} 、{\varGamma _3} 下最可能引起网络结构变化的{\alpha _{ij}} 、{\beta _i} 、{\gamma _j} 排列,避免了无效遍历,减少了求解时间.3. 模拟计算及结果分析
3.1 参数输入
以5个需求站和4个供应站构成的网络进行分析. 供应站列车的到达时刻
{A_{ie}} 、产生的各型空车数量{S_{ ie,\eta }} 见表1. 供应站发出列车f 的最晚编组时刻{L_{if}} 、发往需求站j 和最大空车挂运数量{S_{ if}} 见表2. 需求站发出列车g 的最晚编组时刻{L_{jg}} 及空车需求量{S_{ jg,\chi }} 见表3. 站间空车输送费用{c_{ij}} 、名义旅行时间{t_{ij}} 见表4. 各车站等待车小时费用{c_i} 、{c_j} ,固定作业时间{t_i} 、{t_j} ,空车装车后的发送效益{p_j} 见表5. 规划阶段起始时刻为0,1天为1440 min,表中负号对应前一天的时间. 前一阶段残存空车也按照其实际列车到达情况加入供应站到达列车数据.表 1 供应站列车到达时刻及空车产生数量Table 1. Time of train arrival and number of empty wagons produced in supply stations供应站 产生空车的到达列车编号 1 2 3 4 1 7,5,7,5 412,3,5,5 −439,5,5,5 −100,5,3,5 2 −35,3,5,7 450,3,10,2 −465,7,5,3 −50,5,8,3 3 −367,5,6,5 −497,5,5,5 200,3,4,5 4 −436,5,5,2 −356,5,2,5 80,6,4,5 注:以逗号为分隔的各项依次为列车到达时间/min、平车数目/辆、棚车数目/辆、敞车数目/辆. 表 2 供应站列车最晚编组时刻及最大空车挂运数量Table 2. Train marshalling time limit and maximum number of empty wagons in supply stations供应站 可挂运空车的发出列车编号 1 2 3 4 5 1 222,1,10 627,3,15 −224,5,15 300,2,10 150,4,15 2 203,1,10 693,3,15 −217,5,15 0,2,15 500,4,10 3 −138,2,10 −268,4,10 400,1,10 800,3,10 500,5,10 4 −326,2,15 −246,4,15 200,1,10 400,3,10 340,5,10 注:以逗号为分隔的各项依次为最晚编组时间/min、发往需求站编号、最大空箱挂运数目/辆. 表 3 需求站列车最晚编组时刻及所需空车数量Table 3. Train marshalling time limit and number of empty trains required in demand stations列车
序号需求站编号 1 2 3 4 5 1 660,2,6,2 280,5,2,10 1160,5,8,2 260,2,6,1 400,2,8,2 2 780,5,5,2 560,6,8,5 1400,6,5,2 920,8,5,5 800,10,2,2 注:以逗号为分隔的各项依次为最晚编组时间/min、平车数目/辆、棚车数目/辆、敞车数目/辆. 表 4 站间空车输送费用及名义旅行时间Table 4. Empty wagon transportation cost and nominal travel time between stations供
应
站需求站编号 1 2 3 4 5 1 200,3.3 157,2.3 285,4.6 322,4.8 420,5.9 2 174,3.4 79,1.2 202,3.3 288,4.1 406,5.5 3 298,4.0 172,3.4 186,2.8 193,4.1 204,3.4 4 363,5.5 338,6.2 248,4.9 210,3.5 209,3.1 注:以逗号为分隔的各项依次为空车输送费用/(元•箱)−1、名义旅行时间/h. 3.2 名义模型求解结果
通过MATLAB 2018a调用数学优化软件CPLEX编程求解,得到该供需网络下的最优空车匹配方案,供应站各发出列车挂运的空车来源见表6. 各需求站发出列车的空车来源见表7. 该网络下存在474个第1阶段车流关系,204个第2阶段车流关系,空车输送效益值为53708元.
表 5 车站参数数据Table 5. Station parameter data站点编号 需求站 供应站 pj/(元•箱−1) tj/min cj /(元•h−1) ti/min ci /(元•h−1) 1 673 231 3 205 2 2 485 205 5 225 4 3 598 246 6 218 6 4 659 271 2 105 2 5 805 257 3 表 6 供应站发出列车的空车来源Table 6. Empty wagon source of trains leaving from supply stations供应站编号 列车发出序号 1 2 3 4 5 1 (1,4);(1,6);0 (2,3):0;(2,2) (3,5);(3,5);(3,5) 0;0;0 (4,2);(4,3);(4,1) 2 (1,3);(1,5);(1,2) (2,3);(2,10);(2,2) 0;(3,1);0 (3,7);(3,4);(3,3) 0;0;0 3 0;(1,1);[(1,5),(2,4)] (2,4);(2,5);(2,1) 0;0;0 (1,5);(1,3);0 0;0;0 4 (1,5);(1,5);(1,2) (2,4);(2,2);(2,5) 0;0;(3,2) 0;0;0 [(2,1),(3,6)];(3,2);(3,1) 注:以分号分割的各项依次对应的车型为平车、棚车、敞车; [ ]内为相同车型的空车来源集;( )内各项依次表示到达列车序号、提供空车数量/辆. 以“(1,4);(1,6);0”为例,供应站1发出的第1列列车中,挂运空车的来源为:第1列到达列车提供4辆平车;第1列到达列车提供6辆棚车;该列车中不包含敞车. 表 7 需求站发出列车空车来源Table 7. Empty wagon source of trains leaving from demand station需求站编号 列车发出序号 1 2 1 (2,1,2);(1,1,6);(2,1,2) [(1,1,4),(2,1,1)];(2,1,5);(4,3,2) 2 [(2,4,1),(4,1,4)];(4,1,2);[(2,4,3),(3,1,5),(4,1,2)] (2,4,6);[(2,4,4),(3,1,1),(4,1,3)];[(3,1,4),(4,1,1)] 3 [(1,2,2),(2,2,3)];(2,2,8);(1,2,2) [(1,2,1),(3,4,5)];[(2,2,2),(3,4,3)];(2,2,2) 4 (4,2,2);[(3,2,5),(4,2,1)];(3,2,1) [(1,5,2),(3,2,4),(4,2,2)];[(1,5,4),(4,2,1)];(4,2,5) 5 (1,3,2);[(1,3,7),(2,3,1)];(1,3,2) [(1,3,3),(4,5,7)];(4,5,2);[(1,3,1),(4,5,1)] 注:以分号分割的各项依次对应的车型为平车、棚车、敞车; [ ]内为相同车型的空车来源集合;( )内各项依次表示供应站编号、供应站发出列车序号、空车数量/辆. 以“[(1,1,4),(2,1,1)];(2,1,5);(4,3,2)”为例,需求站1发出的第2列列车中,其装车空箱的来源为:供应站1发出的第1列车提供4辆平车,供应站2发出的第1列车提供1辆平车;供应站2发出的第1列车提供5辆棚车;供应站4发出的第3列车提供2辆敞车. 空车替代情况:需求站2发出的列车2中由列车(4,1,1)供应的1辆敞车为平车用作敞车. 需求站4发出的列车2中由列车(1,5,4)供应的棚车中,有1辆为敞车用作棚车. 需求站5发出的列车1中由列车(1,3,7)供应的棚车中,2辆为敞车用作棚车.
3.3 鲁棒优化模型分析
波动量和波动下限造成的网络结构变化可能导致无可行解,为保证可行解的存在对网络进行下述变化.
1) 为各供应站增加一辆虚拟到达列车,该列车与供应站所有发出列车
f 均存在接续关系,且挂有充足的各类型空车;2) 任一供应站与所有需求站间均增加一列虚拟发出列车,该列车与需求站所有发出列车
g 均存在接续关系且发出列车的空车挂运能力充足;3) 当供应站发出列车空车来源为供应站虚拟到达列车时,供应站库存费用取较大常数,当需求站发出列车车流来源为供应站虚拟发出列车时,需求站发车效益为0.
4) 求得最优解后,除去供应站虚拟到达列车和虚拟发出列车对应车流关系,结合式(1)求解实际效益,作为当前方案下的效益值.
3.3.1 单因素波动影响分析
模型中的不确定因素包括站间旅行时间、供应站及需求站作业时间3项. 各因素的不确定程度由波动量
\hat{t} {{_{ij}}} 、\hat{t} {{_i}} 、\hat{t} {{_j}} 和波动下限{\varGamma _1} 、{\varGamma _2} 、{\varGamma _3} 共同决定,令{q_1}{t_{ij}} = \hat{t} {{_{ij}}} ,{q_2}{t_i} = \hat{t} {{_i}} ,{q_3}{t_j} = \hat{t} {{_j}} . 不同波动下限下效益值随\hat{t} {{_{ij}}} 、\hat{t} {{_i}} 、\hat{t} {{_j}} 的变化即可看作效益值随波动率{q_1} 、{q_{\text{2}}} 、{q_{\text{3}}} 的变化. 单因素波动导致效益值和网络结构的变化如图3所示.由图3可知:
1) 单因素的不确定性导致效益值出现较大变化时,必然伴随车流关系的变化. 只有当减少的车流关系属于当前最优调运方案时效益值才会有较大变化.
2) 相同波动率下效益值差距较大时,波动下限越大效益值越小. 此时,相同波动率下的最优调运方案不同,波动下限越大,网络中存在车流关系越少,网络结构越劣,效益越小.
3) 相同波动率下效益值很接近时,波动下限越大效益值越大. 此时最优调运方案一致,波动下限越大,受到影响的作业或旅行时间越多,车辆的等待时间越小,效益值越大.
4) 从单因素的绝对鲁棒性来看,波动率相同时,效益值越小则受波动率影响越大,故各因素对效益的影响程度为
{t_j} > {t_{ij}} > {t_i} .3.3.2 多因素波动率影响分析
给定
{\varGamma _1}{\text{ = 15}} ,{\varGamma _2}{\text{ = 2}} ,{\varGamma _3} = 3 ,{q_1} 、{q_2} 、{q_3} 之间互动关系如图4.{q_2} 的大小决定第1阶段连接关系的数量,也决定了优化结果的上限. 随着{q_2} 的增大,效益值整体波动幅度逐渐减小,{q_2} 不变时效益值随{q_1} 、{q_3} 的增大而减小.3.3.3 波动下限影响分析
当
q_1 = 5{\text{%}} ,q _2 = 4{\text{%}} ,q_3 = 5{\text{%}} 时,不同{\varGamma _1}、 {\varGamma _2}、{\varGamma _3} 下方案效益值如图5所示.1) 站间旅行时间波动下限影响分析
{\varGamma _2} 、{\varGamma _3} 取值相同,{\varGamma _1} \in [0,2] 时,{\varGamma _1} 对效益有显著影响,并且{\varGamma _1} 越大效益越小.{\varGamma _1} \in [3,20] 时,{\varGamma _1} 对效益无显著影响,{\varGamma _1} 越大效益值越大.2) 供应站作业时间波动下限影响分析
{\varGamma _1} 、{\varGamma _3} 取值相同,{\varGamma _2} 的取值对效益无较大影响,并且{\varGamma _2} 越大效益值越大.3) 需求站作业时间波动下限影响分析
{\varGamma _1} 、{\varGamma _2} 取值相同,{\varGamma _3} \in [0,1] 时的效益明显大于{\varGamma _3} \in [2,5] 时的效益,且效益与{\varGamma _3} 的呈负相关. 但当{\varGamma _3} \in [2,5] 时,虽然效益仍与{\varGamma _3} 呈负相关,但无显著差异.综上,供应站作业时间波动下限对效益值的影响最小. 旅行时间波动下限和需求站作业时间波动下限取值较小时,它们对效益值有较大影响,并且二者间存在相互影响.
4. 结 论
1) 提出供应站到达列车、发出列车以及需求站发出列车间的空车接续时间关系判断方法,并据此分别建立了确定环境与不确定环境下的空车调运优化模型.
2) 考虑构建模型的特征,将确定环境下的模型转化为了线性模型,并结合车流关系的变化是效益值变化的主要原因这一特点,设计了鲁棒优化模型的求解算法,求得了考虑车种替代的空车配流方案.
3) 需求站技术作业时间波动率对效益值影响最大. 供应站技术作业时间波动下限对效益值影响最小. 方案效益值随方案鲁棒性的提高而下降.
-
表 1 模态参数识别结果
Table 1. Modal parameter identification results
序号 方向 模态阶数 频率/Hz 振型 1 横向
弯曲测试结果 4.140 初始 ANSYS
模型模拟结果4.484 2 竖向
弯曲测试结果 4.620 初始 ANSYS
模型模拟结果4.552 3 横向
剪切测试结果 6.895 初始 ANSYS
模型模拟结果7.093 4 纵向
扭转测试结果 8.598 初始 ANSYS
模型模拟结果9.307 5 竖向弯曲 测试结果 10.448 初始 ANSYS
模型模拟结果10.709 表 2 修正参数初始值及取值范围
Table 2. Initial values and range of updated parameters
参数 初始值 下限 上限 数值 θ 值 数值 θ 值 数值 θ 值 E1/Pa 2.00 × 1011 1.00 1.80 × 1011 0.90 2.30 × 1011 1.15 E2/Pa 2.00 × 1011 1.00 1.80 × 1011 0.90 2.30 × 1011 1.15 E3/Pa 2.00 × 1010 1.00 1.75 × 1010 0.88 3.25 × 1010 1.63 ky1/(N·m−1) 1.50 × 107 1.00 4.95 × 106 0.33 2.00 × 107 1.33 kz1/(N·m−1) 1.00 × 108 1.00 5.00 × 107 0.50 1.50 × 108 1.50 ky2/(N·m−1) 1.50 × 107 1.00 4.95 × 106 0.33 2.00 × 107 1.33 kz2/(N·m−1) 1.00 × 108 1.00 5.00 × 107 0.50 1.50 × 108 1.50 ρ2/(kg·m−3) 2.48 × 103 1.00 2.26 × 103 0.91 2.68 × 103 1.08 表 3 ANS修正结果
Table 3. Updated results of ANS method
参数 最大后验概率参数 90% 置信下限 90% 置信上限 数值 θ 值 变化率/% 数值 θ 值 数值 θ 值 E1/Pa 2.19 × 1011 1.10 +10 2.04 × 1011 1.02 2.26 × 1011 1.13 E2/Pa 1.82 × 1011 0.91 −9 1.81 × 1011 0.91 1.89 × 1011 0.95 E3/Pa 3.06 × 1010 1.53 +53 1.94 × 1010 0.97 3.21 × 1010 1.60 ky1/(N·m−1) 0.84 × 107 0.56 −44 0.79 × 107 0.52 1.00 × 107 0.67 kz1/(N·m−1) 0.61 × 108 0.61 −39 0.51 × 108 0.51 1.04 × 108 1.04 ky2/(N·m−1) 1.32 × 107 0.88 −12 1.17 × 107 0.78 1.44 × 107 0.96 kz2/(N·m−1) 1.02 × 108 1.02 +2 0.59 × 108 0.59 1.36 × 108 1.36 ρ2/(kg·m−3) 2.46 × 103 0.99 −1 2.38 × 103 0.96 2.56 × 103 1.03 表 4 NS修正结果
Table 4. Updated results of NS method
参数 最大后验概率参数 90% 置信下限 90% 置信上限 数值 θ 值 变化率/% 数值 θ 值 数值 θ 值 E1/Pa 2.14 × 1011 1.07 +7 2.06 × 1011 1.03 2.14 × 1011 1.14 E2/Pa 1.82 × 1011 0.91 −9 1.82 × 1011 0.91 1.97 × 1011 0.99 E3/Pa 3.20 × 1010 1.60 +60 2.19 × 1010 1.10 3.20 × 1010 1.61 ky1/(N·m−1) 0.83 × 107 0.55 −45 0.81 × 107 0.54 0.97 × 107 0.65 kz1/(N·m−1) 0.62 × 108 0.62 −38 0.53 × 108 0.53 1.03 × 108 1.03 ky2/(N·m−1) 1.32 × 107 0.88 −12 1.22 × 107 0.81 1.53 × 107 1.02 kz2/(N·m−1) 1.06 × 108 1.06 +6 0.61 × 108 0.61 1.38 × 108 1.38 ρ2/(kg·m−3) 2.42 × 103 0.98 −2 2.38 × 103 0.96 2.60 × 103 1.05 表 5 修正后模态频率对比
Table 5. Comparison of updated modal frequencies
模态频率/Hz 测试值/Hz 初始值/Hz 初始误差/% NS 修正值/Hz NS 误差/% ANS 修正值/Hz ANS 误差/% 1 4.140 4.484 8.30 4.162 0.53 4.159 0.47 2 4.620 4.552 −1.47 4.627 0.15 4.617 −0.06 3 6.895 7.093 2.87 6.930 0.51 6.938 0.63 4 8.598 9.307 8.24 8.580 −0.21 8.552 −0.54 5 10.448 10.709 2.50 10.492 0.42 10.463 0.15 注:相对误差=(修正值−测试值)/测试值. -
[1] MEI Z, WU B, BURSI O S, et al. Hybrid simulation with online model updating: application to a reinforced concrete bridge endowed with tall piers[J]. Mechanical Systems and Signal Processing, 2019, 123: 533-553. doi: 10.1016/j.ymssp.2019.01.009 [2] CHEN Z, SUN H. Sparse representation for damage identification of structural systems[J]. Structural Health Monitoring, 2020, 20(4): 1644-1656. [3] MO J, WANG L, GU K X. A two-step interval structural damage identification approach based on model updating and set-membership technique[J]. Measurement, 2021, 182: 109464.1-109464.19. [4] PU Q H, HONG Y, CHEN L J, et al. Model updating-based damage detection of a concrete beam utilizing experimental damped frequency response functions[J]. Advances in Structural Engineering, 2019, 22(4): 935-947. doi: 10.1177/1369433218789556 [5] ZENG J C, KIM Y H. Identification of structural stiffness and mass using bayesian model updating approach with known added mass: numerical investigation[J]. International Journal of Structural Stability and Dynamics, 2020, 20(11): 2050123.1-2050123.31. [6] ASTROZA R, ALESSANDRI A, CONTE J P. A dual adaptive filtering approach for nonlinear finite element model updating accounting for modeling uncertainty[J]. Mechanical Systems and Signal Processing, 2019, 115: 782-800. doi: 10.1016/j.ymssp.2018.06.014 [7] 万华平,任伟新,黄天立. 基于贝叶斯推理的随机模型修正方法[J]. 中国公路学报,2016,29(4): 67-76,95. doi: 10.3969/j.issn.1001-7372.2016.04.009WAN Huaping, REN Weixin, HUANG Tianli. Stochastic model updating approach by using Bayesian inference[J]. China Journal of Highway and Transport, 2016, 29(4): 67-76,95. doi: 10.3969/j.issn.1001-7372.2016.04.009 [8] JOUBERT D J, MARWALA T. Monte Carlo dynamically weighted importance sampling for finite element model updating[M]//Conference Proceedings of the Society for Experimental Mechanics Series. Cham: Springer International Publishing,2016: 303-312. [9] 彭珍瑞,郑捷,白钰,等. 一种基于改进MCMC算法的模型修正方法[J]. 振动与冲击,2020,39(4): 236-245.PENG Zhenrui, ZHENG Jie, BAI Yu, et al. A model updating method based on an improved MCMC algorithm[J]. Journal of Vibration and Shock, 2020, 39(4): 236-245. [10] 秦世强,廖思鹏,黄春雷,等. 基于自适应Kriging模型的人行斜拉桥有限元模型修正[J]. 中山大学学报(自然科学版)(中英文),2021,60(6): 43-53.QIN Shiqiang, LIAO Sipeng, HUANG Chunlei, et al. Adaptive Kriging model based finite element model updating of a cable-stayed pedestrian bridge[J]. Acta Scientiarum Naturalium Universitatis Sunyatseni, 2021, 60(6): 43-53. [11] 王其昂,王腾,倪一清. 基于卷积神经网络的高铁车轮损伤识别方法研究[J]. 中国矿业大学学报,2020,49(4):781-787.WANG Qi’ang, WANG Teng, NI Yiqing. Damage detection of wheels for high-speed rail based on convolutional neural network[J]. Journal of China University of Mining & Technology, 2020, 49(4):781-787. [12] SKILLING J. Nested sampling for general Bayesian computation [J]. Bayesian Analysis,2006,1:833-859. [13] HANDLEY W J, HOBSON M P, LASENBY A N. Polychord: next-generation nested sampling[J]. Monthly Notices of the Royal Astronomical Society, 2015, 453(4):4384-4398. [14] CAO T T, ZENG X K, WU J C, et al. Integrating MT-DREAMzs and nested sampling algorithms to estimate marginal likelihood and comparison with several other methods[J]. Journal of Hydrology, 2018, 563: 750-765. doi: 10.1016/j.jhydrol.2018.06.055 [15] HIGSON E, HANDLEY W J, HOBSON M P, et al. Dynamic nested sampling: an improved algorithm for parameter estimation and evidence calculation[J]. Statistics and Computing, 2019, 29(5): 891-913. doi: 10.1007/s11222-018-9844-0 [16] TRASSINELLI M, CICCODICOLA P. Mean shift cluster recognition method implementation in the nested sampling algorithm[J]. Entropy, 2020, 22(2): 185. doi: 10.3390/e22020185 [17] QIAN F, ZHENG W. An evolutionary nested sampling algorithm for Bayesian model updating and model selection using modal measurement[J]. Engineering Structures, 2017, 140: 298-307. doi: 10.1016/j.engstruct.2017.02.048 [18] XU X K, HONG Y, CHEN L J, et al. Hybrid nested sampling method for identifying the uncertainty of the high-dimensional updating parameters in Bayesian structural model updating[J]. Advances in Structural Engineering, 2022, 25(8): 1730-1744. doi: 10.1177/13694332211069511 [19] 王坤阳,公茂盛,左占宣. 基于贝叶斯理论嵌套抽样的结构物理参数识别研究[J]. 振动与冲击,2022,41(7): 74-80.WANG Kunyang, GONG Maosheng, ZUO Zhanxuan. Structural physical parameter identification based on Bayesian theory and nested sampling[J]. Journal of Vibration and Shock, 2022, 41(7): 74-80. [20] 杨朋超,薛松涛,谢丽宇. 消能减震建筑结构的贝叶斯有限元模型修正[J]. 土木工程学报,2021(S01): 13-19,47.YANG Pengchao, XUE Songtao, XIE Liyu. Bayesian finite-element model updating of passively controlled building structures[J]. China Civil Engineering Journal, 2021(S01): 13-19,47. [21] NEAL R M. Slice sampling[J]. The Annals of Statistics, 2003, 31(3): 705-767. [22] DONG X J, WANG Y. A comparative study of frequency-domain finite element updating approaches using different optimization procedures[C]//8th European Workshop on Structural Health Monitoring (EWSHM 2016). Bilbao: [s.n.],2016: 1-10. [23] 洪彧. 基于振动信号的桥梁结构模态参数识别与模型修正研究[D]. 成都:西南交通大学,2019. -