基于响应面法的UHPC-NC连续箱梁桥优化设计.pdf

762023年6 月中路夕第4 3 卷第3 期DOI:10.14048/j.issn.1671-2579.2023.03.011基于响应面法的UHPC-NC连续箱梁桥优化设计关维溢，段璐琦（1.佛山市建盈发展有限公司，广东佛山528315；2.长沙理工大学土木工程学院，湖南长沙410114)摘要：为了解决大跨桥梁后期出现的跨中下挠和开裂两大问题，提出UHPC-NC混合梁设计方案，并基于响应面法对该设计部分设计参数进行优化。通过对局部参数即跨中截面的顶板、腹板和底板厚度和总体参数如UHPC替换比例、梁体曲线次数和跨中梁高等参数进行优化，最终获得最佳参数组合为顶板厚度2 2 cm、腹板厚度2 5cm、底板厚度3 2cm、U H PC替换比例0.6 3、曲线次数1.8、跨中梁高3.3 5m。然后以最佳参数组合为基础建立有限元模型与原设计模型进行比较，结果表明各墩墩顶弯矩减小了12%以上，全桥弯曲应变能下降了3 4.7%，各关键截面的上下缘应力普遍减小了12%，最大挠度降低了2 3.3%。关键词：大跨连续梁桥；UHPC；响应面法；结构优化中图分类号：U443.35文献标志码：A0引言对现存的连续梁桥跟踪调研发现，此类桥梁在长期运营过程中常出现跨中下挠、梁体开裂两大病害，而受限于普通混凝土材料性能以及连续梁桥本身固有结构体系，这两大病害并不能得到很好的解决1。为解决这两大病害，提出了UHPC-NC连续混合箱梁的新方案2 。即对跨中区域的普通混凝土梁体采用UHPC梁体进行替换，同时为了充分发挥UHPC的材料性能，对桥梁的部分结构参数进行优化以获得最优化参数组合。但是目前对于结构参数的优化往往都是针对各个设计变量进行单独分析，无法考虑到各个设计变量之间的相互影响。贾丽君等3 采用控制变量法对单一设计参数进行改变，实现了UHPC连续刚构桥的设计优化，并取得了良好的效果；宋恒扬等4 对一座空腹式连续刚构桥的部分重要设计参数如梁底曲线次数和下弦梁高等，采用正交试验设计后进行了优化分析，获得了最佳经济性。以上的优化设计方法，常常没有考虑到各个因素之间的耦合影响，导致优化效果不够理想，且优化结果的精度也不够高。响应面法(RSM)在生化领域5-7 中有着广泛的应用，相较于其他优化方法，该法能直观显示多个变量与目标之间的相关性，且收敛性更加优越，但这种方法需要严格选取合适的响应值以保证优化结果的有效性及准确性，这就需要一定的工程设计经验。目前越来越多的工程师将响应面法应用于土木工程中8-9 ,本文将其引人到桥梁结构优化设计中以解决上述病害。1基基于响应面法的设计参数优化方法基于响应面的优化理论其主要内容包括试验设计、参数筛选、响应面模型建立、精度检验以及求解最优化参数。土木工程中常见的试验设计方法有：中心复合试验设计（CCD）、Bo x-Be h n k e n 设计（BBD）和D-optimal设计10 等。本文采用BBD试验设计方法，该试验设计是一种不包含嵌入影响因素或者部分影响因素的试验设计方法，常用于研究影响因素与响收稿日期：2 0 2 2-0 5-10（修改稿）基金项目：国家重点研发计划项目（编号：2 0 19 YFC1511002）；国家自然科学基金资助项目（编号：517 7 8 0 6 8，52 0 7 8 0 58）；湖南省自然科学基金创新研究群体项目（编号：2 0 2 0 JJ1006）；湖南省教育厅自然科学研究重点项目（编号：2 1A0196）；长沙市自然科学基金资助项目（编号：kq2202209）作者简介：关维溢，男，大学本科，工程师.E-mail:772023年第3 期关维溢2面法的UHP梁桥优化设计人应结果的非线性关系且存在至少3 个影响因素的试验。在该试验设计中，选取N维超立方体每个边的中点作为抽样点，每个影响参数设置了一1、0、1三个水平，分别代表最小值、中心点以及最大值。在进行响应面模型拟合前，应该对所选择的参数因子进行筛选，确定其中显著的参数因子，淘汰不显著的参数因子，避免后续试验中由于参数太多，使得计算过于复杂，浪费试验资源。目前在响应面方法中，常采用方差分析法即F检验法进行参数显著性检验。F检验法表达式：FA=SSA/fAF(fA.f.)(1)SS./f.式中：SSA为各因素引起的偏差平方和；SS.为试验误差引起的偏差平方和；fA为各因素的自由度；f.为偏差的自由度。假定显著水平临界值取0.0 1,当FFi-0.01(Ja,f),则认为该因素影响非常显著；当F1-0.05(fa,f.)FFi-0.0(fA,f.),则认为该因素影响显著；当FF1-0.1(fA,f.),则认为该因素影响不显著。常用的响应面函数有多种类型，如：多项式函数、幂函数、径向基函数、BP神经网络等12 ,土木工程中一般采用多项式的函数形式，该函数在满足精度要求的情况下同时具备良好的可操作性。常用的多项式函数的形式有一阶多项式响应面函数、二阶多项式响应面函数和三阶多项式响应面函数10 。由于一阶多项式响应面函数只能描述线性关系，无法描述实际工程中较为复杂的响应关系，而三阶多项式响应面函数项及待定系数数量过多不便于实际使用，故本文选择二阶多项式响应面函数。二阶线性响应面函数表达式：a.+2之7y=+iai(2)2i=1i=1j=1i=1式中：y为响应值;i,为待优化参数值;oviv和it为响应面函数的待定系数，这些系数通过样本试验确定。在取得最终优化结果前，需要对响应面函数的精度进行检验，以保证结果的有效性13 。本文采用R?准则对拟合的响应面函数精度进行检验。Zy.(i)-y(i)R2=1-豆(3)Z(i)-习=1式中：m为试验样本点的数量；y（i)为第i个试验样本点的响应面模型计算值；(i）为试验样本点的有限元模型计算值；为y(i)平均值。如果R的值均趋近于1,则表明响应面函数的拟合情况良好。2UHPC-NC混合梁桥设计参数优化2.1工程概况以某特大桥为工程背景，桥上部构造为五跨预应力混凝土连续箱梁，全桥主梁采用C55混凝土，三向预应力体系，桥梁立面图如图1所示。梁底曲线次数为2，箱梁根部梁高6.9 m，跨中梁高3.0 m，顶板厚28cm，跨中腹板厚4 0 cm，底板厚3 2 cm。箱梁顶板的横向宽度为12.7 5m，箱底宽度为7.0 m，翼缘悬臂长度为2.8 7 5m，关键位置横截面图如图2 所示。4640006.7001100011000110006700口口口图1某特大桥立面布置图（单位：cm）1275/21275/20084006970支点截面跨中截面图2支点及跨中横截面图（单位：cm）2.2优化目标由于桥梁的设计包含众多参数，本文受限于篇幅，只对部分重要的局部参数以及整体参数分别进行优化设计。考虑是对跨中部分区段进行UHPC材料替换，并进行结构设计优化，故局部参数首选对跨中截面的尺寸一一顶板厚度、底板厚度和腹板厚度进行优化，以达到最好的减重效果。然后对总体参数而言，选取UHPC的替换长度、跨中梁高和梁底曲线次数进行试优化2.3局部参数优化试验2.3.1跨中截面尺寸参数选择将顶板厚度、腹板厚度以及底板厚度3 个设计变78公中路外第4 3 卷量分别用A、B和C代替。依据参考文献13 并结合工程实际确定各个设计变量范围如表1所示，表1各设计变量范围取值范围A/cmB/cmC/cm最小值202525最大值2840322.3.2优化评价标准选择在局部参数优化中主要考虑两个评价参数作为响应值，分别为：(1)Ap/(EI)在自重荷载作用下，主梁的最大挠度为：gL4.Ap(4)maxkEIEI式中：g为重力加速度；l为主梁跨度；A为主梁截面面积；为主梁材料密度；k为主梁边界条件与跨径布置决定的系数；EI为截面抗弯刚度。在桥梁材料和结构形式已经确定时，A/（EI)反映了结构刚度受自重荷载影响的大小，其值越低表示结构主梁刚度越高。在桥梁材料和结构形式已经确定时，Ap/EI)反映了结构刚度受自重荷载影响的大小，其值越低表示结构主梁刚度越高。（2）截面抗弯效率参数n6=(K,+K,)/h(5)式中：K为截面上核心距K，为截面下核心距；h为梁截面高度；n反映截面对荷载的综合抗弯能力，仅与截面形式有关，其值越大表明截面抗弯效率越高。取Ri=Ap/EI)和R2=n这两个评价参数作为响应值，从而选出最优化设计参数组合2.3.3响应面模型的建立与分析在确定好各个待优化参数的选择范围后，采用BBD试验设计方案，进行三因素三水平的响应面试验，中心点重复3 次，共进行15次试验，并计算得出各个样本点评价指标结果如表2 所示。表2试验设计方案与结果顶板厚度腹板厚度底板厚度组数R1R2/102(A)/cm(B)/cm(C)/cm12440.025.052.323760.333 622432.528.549.973962.573 832032.532.048.357163.385442432.528.549.973962.573852025.028.548.070364.806 7续表2顶板厚度腹板厚度底板厚度组数RiR2/102(A)/cm(B)/cm(C)/cm62032.525.050.458462.424472432.528.549.973962.573882832.532.049.745162.426792832.525.052.204261.3452102425.025.050.268 863.6892112425.032.047.912.064.6841122440.032.050.086961.3670132840.028.551.804.960.4980142825.028.549.789963.6221152040.028.550.371 461.3276图3 为各参数对响应值的显著性图。图3 中实线表示P=0.0001。由图3 可知：单个因素A、B与C和交互相AB和AC以及二次项C2对响应值Ri的影响极为显著，交互相BC和二次项A?的影响较为显著，二次项B2对响应值的影响并不显著；单个因素A、B和C和交互相AB以及二次项A、B2 和C2对响应值R2的影响极为显著，交互相AC的影响较为显著，交互相BC对响应值的影响并不显著。0.7R0.6一R20.50.4P0.30.20.10一ABCABACBCA2B22图3各个参数对响应值的显著性图对试验结果进行二次多元回归拟合，并舍去部分影响较小的参数，得到Ri和R2关于顶板厚度A、腹板厚度B和底板厚度C3个因素的多元二次回归模型方程分别为：Ri=57.919+0.338A+0.172B-1.039C0.002AB-0.006AC+0.001BC+0.002A2+0.015C2R2=(64.033-0.084A-0.368B+0.602C+0.003AB+0.002AC-0.004A2+792023年第3 期关维溢，面法的UHPC桥优化设计0.001B2-0.009C2)/100两个模型的校正决定系数R均为1,表明模型与实际情况拟合很好。针对回归模型，运用Design-Expert12软件中的Optimization功能，得到跨中截面的最优化设计参数，顶板厚度2 2.0 2 cm、腹板厚度2 5.3 4 cm、底板厚度3 1.7 cm为最优参数。为了便于在实际工程中的应用以及标准化设计，将最优化设计参数分别修正为2 2 cm、2 5c m 和3 2 cm。2.4总体参数优化试验2.4.1参数比选及范围在优化了局部参数的基础上，为了充分发挥UHPC材料的性能，还需要对部分总体设计参数进行优化。选取UHPC的替换长度与主跨跨径的比值k、梁底曲线次数和跨中梁高作为待优化目标，并用字母D、E和F代替表示。通过查阅相关文献13 并结合工程实际确定各个待优化设计参数变化范围如表3所示。表3各设计参数变化范围取值范围k(D)梁底曲线次数（E）跨中梁高（F)/m最小值0.11.32最大值0.92.542.4.2优化评价标准选择在确定了设计参数的变化范围后，还需要选取合适的评价标准作为响应值，从而得出最优化设计，此次优化共选取了5个参数作为评价标准。(1）弯曲应变能R可以用主梁的弯曲应变能来评价主梁弯矩分布的均匀性，其值越小，表明整个主梁的弯矩分布越均匀，主梁受力状态越好。对离散的杆系结构，该式可以写成：LU.-27(Mi+M)(6)4(EI)i=1式中：m为主梁单元数；L和（EI)，分别为主梁第i单元的长度和抗弯刚度，对于变截面连续梁桥，单元两侧的截面惯性矩不同，采用前文的等效截面惯性矩进行计算；MiL和Mi分别为主梁第i单元的左端弯矩和右端弯矩。（2）跨中挠度R2以恒载作用下的跨中最大挠度作为评价指标。跨中挠度越小越好(3）U H PC 用量R3对于采用UHPC材料的连续梁桥，UHPC不仅本身材料价格大大高出普通混凝土，而且其施工难度也较普通混凝土更大，UHPC的使用方量会直接影响工程造价，故可以在取得最好的优化结构的同时尽量减少UHPC材料的使用量。（4）跨中最大弯矩R4以自重作用下跨中截面的最大弯矩作为评价标准，跨中弯矩越小越好（5）支点最大弯矩R5以自重作用下支点截面的最大弯矩作为评价标准，支点弯矩越小越好2.4.3响应面模型的建立与分析在确定好各个待优化参数的选择范围后，进行三因素三水平的响应面试验，中心点重复3 次，共进行15次试验。试验设计方案及计算结果如表4 所示。图4 为各参数对响应值的显著性图。表4试验设计方案与结果曲线次跨中梁高弯曲应变能跨中度UHPC用量跨中最大弯矩支点最大弯矩组数(D)数(E)(F)/mR./kJR2/mmR:/m3R4/(kNm)Rs/(kNm)10.51.341215.5465.6053662.593598.6-36108920.91.921636.93133.8788757.249411.6-31001530.51.931402.6089.1550175.870669.0-34836240.12.531.775.95105.609982.877220.3-36640150.11.941585.5877.7110608.2102952.0-39038960.11.921966.87135.369370.449551.2-37021570.52.521834.13144.7446789.850.337.1-33246680.51.931402.6089.1550175.870669.0-348.36280公中路外第4 3 卷续表4曲线次跨中梁高弯曲应变能跨中挠度UHPC用量跨中最大弯矩支点最大弯矩组数(D)数(E)(F)/mRi/kJR2/mmR3/m3R4/(kNm)Rs/(kNm)90.51.321529.51106.0848646.247601.2-360941100.11.331587.1090.3010019.279.542.3-401544110.51.931402.6089.1550175.870669.0-348 362120.92.531415.15102.0991280.174514.3-303286130.91.941161.3273.4496330.794954.9308536140.91.331259.2286.0894479.974340.6-316199150.52.541272.2772.8152563.893 414.3-343289图4 中实线表示P=0.0001。由图4 可知：单次相D、E和F和二次项D、F2 对响应值R的影响极为显著，交互相DF、EF的影响较为显著，而交互相DE和二次项E2对Ri的影响并不显著；单次相F对响应值R2的影响极为显著，单次相E、交互相EF和二次项D、F2 的影响较为显著，而单次相D、交互相DE、DF和二次项E2对R2的影响并不显著；单次相D、F对响应值R3的影响极为显著，单次相E、交互相DF、DF、二次项D?的影响较为显著，而交互相EF和二次项E2、F2 对R3的影响并不显著；单次相F对响应值R4的影响极为显著，单次相D和二次项D?的影响较为显著,其余参数对R4的影响并不显著；单次相D、E对响应值Rs的影响极为显著，单次相F、交互相DE、DF、EF和二次项D?的影响较为显著，而二次项E?和F2对Rs的影响并不显著。1.0R0.9R20.8R30.7R40.6Rs0.50.40.30.20.10-0.1DEFDEDFEFD2E2F2图4各个参数对响应值的显著性图对试验结果进行二次多元回归拟合，并舍去部分影响较小的参数，得到各个响应值关于参数UHPF替换比例D、曲线次数E和跨中梁高F3个因素的多元二次回归模型方程分别为：Ri=2002.58-932.527D+569.413E-406.628F-58.989DF-103.29EF+723.654D2+69.295F2R2=154.938+59.877E54.75F-13.104EF+45.797D+8.618F2R:=2209.993+91165.951D2772.494E327.183F-3295.275DE+3959.835DF+6 611.581D2R4=4020.35322454.563D+33124.062F-4911.206DF+27234.703D2Rs=-4.536X105+7.493X10*D+6.121X10*E-4.977X103F-2.316X10*DE+1.35X10*DF-4.448X10EF+1.926X10D?5个模型的校正决定系数R分别为0.9 9 9 2、0.9962、1、0.9 9 6 7、0.9 9 8 7，表明模型与实际情况拟合很好。针对回归模型，运用Design-Expert12软件中的Optimization功能，得到总体参数的最优化设计组合为：k=0.6317、曲线次数为1.7 6 4、跨中梁高为3.346m。为了便于在实际工程中的应用以及标准化设计，将最优化设计参数分别修正为k一0.6 3、曲线次数取1.8、跨中梁高为3.3 5m。3优化结果对比3.1主梁弯矩对比全桥关键截面弯矩如表5所示。由表5可以看出：相比于原设计，优化设计后的模型中跨跨中位置弯矩有小幅度增大，由53 6 9 9kNm增大到了550 2 0 kNm，弯矩增大了2.5%。812023年第3 期关维溢，应面法的UHP梁桥优化设计车续租表5优化前后各关键截面弯矩弯矩/(kNm）截面位置优化前优化后边跨1/43297533093边跨1/25865089边跨3/4-119548-96888边墩墩顶337428-293117次中跨1/4-9005134次中跨1/27712977217次中跨3/4-32435-22062中墩墩顶384 288-337510中跨1/4-44150-33160中跨1/25369955020而位于两个墩顶位置处的负弯矩有大幅减小，其中边墩墩顶位置的负弯矩由一3 3 7 4 2 8 kNm减小到了-2 9 3 117 kNm，减小了13.1%，中墩墩顶位置的弯矩由-3 8 4 2 8 8 kNm减小到了-3 3 7 510 kNm，减小了12.2%。虽然跨中区段的弯矩有了小幅度增加，但考虑到该区域采用的是高抗拉抗压的UHPC材料，故小幅度的弯矩增幅对该局部区域的影响有限。而墩顶区域的弯矩大幅度减小则能很大程度上改善该区域梁段受力，对于还是使用普通混凝土材料的墩顶区域，降低该区域的内力水平能有效减小其徐变程度，从而进一步改善跨中挠度。同时通过计算得出，在自重作用下主梁的弯矩应变能由原设计的17 8 7 kJ下降到了116 7 kJ，下降了34.7%，由此可见，优化后的设计比原设计的弯矩分布更加均匀。3.2主梁应力对比全桥关键截面上下缘应力如表6 所示。表6优化前后各关键截面应力优化前/MPa优化后/MPa截面位置上缘下缘上缘下缘边跨1/4-3.024.84-2.854.53边跨1/21.02-1.480.58-0.95边跨3/45.967.224.625.54边墩墩顶7.998.776.89-7.55次中跨1/40.710.96-1.141.76次中跨1/27.1311.50-6.7110.70续表6优化前/MPa优化后/MPa截面位置上缘下缘上缘下缘次中跨3/42.13-2.991.532.47中墩墩顶9.22-10.108.05-8.81中跨1/41.51-2.180.95-1.54中跨1/2-4.998.06-4.807.64由表6 可以得出：优化后模型关键截面的上下缘应力均有不同程度的降低。边墩墩顶的上缘应力由7.99MPa降低为6.8 9 MPa，降低了13.8%；下缘应力由-8.7 7 MPa降低到一7.55MPa，降低了13.9%。次中跨跨中的上缘应力由一7.13 MPa降低为一6.7 1MPa，降低了5.9%；下缘应力由11.5MPa降低为10.7MPa，降低了7%。中墩墩顶的上缘应力由9.2 2MPa降低为8.0 5MPa，降低了12.7%；下缘应力由一10.1MPa降低为一8.8 1MPa，降低了12.8%。中跨跨中位置的上缘应力由一4.9 9 MPa降低到一4.8MPa，降低了3.8%；下缘应力由8.0 6 MPa降低到7.6 4MPa，降低了5.2%。可以发现虽然跨中处的应力变幅幅度较小，但墩顶处梁单元的上下缘应力都有了大幅降低。3.3主梁挠度对比全桥关键截面挠度如表7 所示。由表7 可得出：全桥最大度出现于次中跨跨中位置。原设计最大挠度为一10 3 mm，优化后模型最大挠度为一7 9 mm，降低了2 3.3%，说明采用UHPC材料替换跨中部分梁段并进行优化设计后能有效降低跨中挠度。表7优化前后各关键截面挠度挠度/mm截面位置优化前优化后边跨1/4-9.7-8.1边跨1/2一2.13.0边跨3/47.45.1次中跨1/4-56.144.1次中跨1/2-103.0-79.0次中跨3/4一4 8.0-37.9中跨1/4-15.214.4中跨1/2-48.4一4 0.582公中路第4 3 卷外4结论为了解决连续梁桥在长期运营过程中跨中过度下挠，跨中梁段开裂的问题，本文提出一种采用UHPC材料替换跨中部分普通混凝土梁段的设计，并基于响应面法对该设计方案的部分设计参数进行优化，在获得最好经济效益的基础上充分发挥UHPC材料的力学性能。得到如下结论：（1）在局部参数优化中，选取跨中截面的顶板厚度、腹板厚度和底板厚度作为待优化参数，并选取Ap/（EI)和截面抗弯效率参数n两个参数作为优化评价指标，通过响应面分析得出最优化设计参数为顶板厚度为2 2 cm、腹板厚度为2 5cm和底板厚度为32cm。（2）在总体参数优化中，选取UHPC替换比例k、梁高变化曲线次数和跨中梁高作为待优化参数，选取弯曲应变能、挠度、UHPC用量和跨中、支点最大弯矩作为优化评价指标，通过响应面分析得到最优化设计参数为k一0.6 3、曲线次数取1.8、跨中梁高为3.35m。（3）采用最优化参数建立有限元模型，并选取截面弯矩、应力和挠度等参数作为指标与原设计方案进行比较。结果显示：采用最优化设计的墩顶位置弯矩较原设计有大幅减小，基本减小12%以上，同时弯曲应变能降低了3 4.7%，表明最优化设计的弯矩分布比原设计更加合理。同时梁单元的上下缘应力也有不同幅度的减小，特别是墩顶截面的应力减小12%以上，这对于还是采用普通混凝土材料的墩顶梁段受力更加良好。通过比较最优化设计和原设计两个模型的全桥挠度可以发现，采用新方案的最大挠度降低了2 3.3%，表明优化后的方案在降低跨中挠度这一方面具有非常好的效果。参考文献：1 李庆海.高速公路混凝土梁桥常见病害总结及成因分析J.上海公路,2 0 2 0(4):6 0-6 4,7 0.2 赵佳兴.大跨预应力UHPC-NC连续混合箱梁桥分析与优化设计D.长沙：湖南大学,2 0 18.3 贾丽君，林赞笔，袁勇根，等.预应力UHPC连续刚构桥的优化设计J.沈阳工业大学学报,2 0 17,3 9(5):59 1-59 5.4宋恒扬，胡俊，彭元诚.空腹式连续刚构桥设计参数的正交试验研究J.世界桥梁,2 0 15,4 3(1):55-58.5 权帆，朱文秀，张晴，等.响应面法优化藜麦发酵酸奶的工艺研究J.食品研究与开发，2 0 2 2,4 3(8):13 3-13 9.6贾磊，程海鹰，魏玮，等.响应面法的泡沫沥青发泡控制参数优化J.机械设计与制造，2 0 2 2,3 7 2(2)：7 1-7 3,7 7.7FERDOSIAN I,CAMOES A.Eco-efficient ultra-highperformance concrete development by means of responsesurface methodologylJ.Cement and Concrete Composites,2017,84:146-156.8DEY T K,MUKHOPADHYAY T,CHAKRABARTI A,et al.Efficient lightweight design of FRP bridge deckJ.Proceedings of the Institution of Civil Engineers-Structuresand Buildings,2015,168(10):697-707.9XU L,BI K M,GAO J F,et al.Analysis on parameteroptimizationofdampersSoflong-spandouble-towercable-stayedbridgesJ.StructureandInfrastructureEngineering,2020,16(9):1286-1301.10徐智杰,郭细伟.双肢薄壁墩连续刚构桥的参数优化J.工程与建设,2 0 2 0,3 4(5):8 50-8 52,8 58.11刘晓慧.基于索力响应面的斜拉桥多尺度有限元模型修正D.石家庄：石家庄铁道大学,2 0 17.12秦仙蓉，潘杰，徐俭，等.塔式起重机结构有限元模型修正的响应面方法J.振动与冲击,2 0 18,3 7(6):2 4 4-2 50.13 梁乾敏.连续刚构桥的合龙顶推模型修正与顶推力线性规划解法研究D.西安：长安大学,2 0 2 1.