数学建模教材(第四章)

第4章数学规划模型本章研究数学规划模型，其中包括：线性规划、整数规划、非线性规划、多目标规划与动态规划等内容．线性规划模型线性规划是运筹学的一个重要分支，随着计算机技术的发展，线性规划不仅在理论上已趋向成熟，而且在实际应用中也日益广泛与深入．本节将借助Lingo数学软件对线性规划模型进行求解．4.1.1问题的提出在生产管理和经营活动中经常提出一类问题，即如何合理地利用有限的人力、物力、财力等资源，以便得到最好的经济效果．引例1 普通生产计划安排问题某工厂在计划期内要安排生产Ⅰ、Ⅱ两种产品，已知生产单位产品所需的设备台时及A、B两种原材料的消耗，如表4-1所示．该工厂每生产一件产品Ⅰ可获利2元，每生产一件产品Ⅱ可获利3元，问应该如何安排计划使该工厂获利最多表普通生产计划安排问题ⅠⅡ设备原材料A 原材料B利润1422438台时16kg12kg引例2 奶制品的生产计划问题一奶品加工厂用牛奶生产A、B两种奶制品，1桶牛奶可以在甲类设备上用12小时加工成3公斤A，或者在乙类设备上用8小时加工成4公斤B，根据市场需求，生产的A、B全部能售出，且每公斤A获利24元，每公斤B获利16元．现在加工厂每天能得到50桶牛奶的供应，每天正式工人总的劳动时间为480小时，并且甲类设备每天最多能加工100公斤A，乙类设备的加工能力没有限制．试为该厂制定一个生产计划，使每天获利最大，并进一步讨论以下３个附加问题：⑴若用35元可以买到1桶牛奶，应否做这项投资若投资，每天最多购买多少桶牛奶⑵若可以聘用临时工人以增加劳动时间，付给临时工人的工资最多是每小时几元⑶由于市场需求变化，每公斤A的获利增加到30元，应否改变生产计划4.1.2模型建立1.引例1普通生产计划安排问题的模型建立对于引例1，可以设x、y分别表示在计划期内产品Ⅰ、Ⅱ的产量．若用z表示利润，这时，23z x y =+．因为设备的有效台时为8，因此应有限制条件：28x y +≤；同理考虑原材料的不同限制，可得如下限制条件：416x ≤，412y ≤．综上所述，该计划问题可用数学模型表示为：目标函数：max 23z x y =+约束条件：28416412x y x y +≤⎧⎪≤⎨⎪≤⎩0,0x y ≥≥2.引例2奶制品生产计划问题的模型建立对于引例2，可以设每天用1x 桶牛奶生产A ，用2x 桶牛奶生产B ．类似引例1可得奶制品生产计划问题的数学模型：目标函数：12max 7264z x x =+约束条件：12121501284803100x x x x x +≤⎧⎪+≤⎨⎪≤⎩120,0x x ≥≥从以上两例可以看出，他们都属于同一类优化问题，他们的共同特征： ⑴每一个问题都用一组决策变量12(,,,)n x x x 表示某一方案；这组决策变量的值就代表一个具体方案，一般这些变量取值都是非负的；⑵存在一定的约束条件，这些约束条件可以用一组线性等式或线性不等式来表示；⑶都有一个要求达到的目标，它可用决策变量的线性函数来表示，这个函数称为目标函数．满足以上三个条件的数学模型称为线性规划数学模型．其一般形式为：目标函数：1122max(min)n n z c x c x c x =+++约束条件：11112211211222221122(,)(,)(,)n n n n m m mn n ma x a x a xb a x a x a x b a x a x a x b ++≤=≥⎧⎪++≤=≥⎪⎨⎪⎪++≤=≥⎩120,0,0n x x x ≥≥≥4.1.3模型求解1.引例1普通生产计划安排问题的模型求解使用Lingo 数学软件对引例１进行求解，编写程序如下：max =2*x+3*y; x+2*y<8; 4*x<16; 4*y<12; end求解结果为：目标函数的最大值14z =，即可获得最大利润１４元；最优生产计划方案是：生产产品Ⅰ4件，生产产品Ⅱ2件．2.引例2奶制品生产计划问题的模型求解使用Lingo 数学软件对引例2进行求解，编写程序如下：max =72*x1+64*x2;x1+x2<50;12*x1+8*x2<480; 3*x1<100; end求解结果为：每天用20桶牛奶生产A ，用30桶牛奶生产B ，最大利润是3360元．下面来回答三个附加问题：⑴针对附加问题1，可假设应投资购买x 桶牛奶，目标函数应修改为：12max 726435*z x x x =+-关于牛奶的约束条件也应作相应的修改：1250x x x +<+通过编程求解得：最大利润增加到3490元，因此，应作该项投资，每天最多购买10桶牛奶．⑵针对附加问题2，首先将劳动时间480小时增加1个小时，对原问题进行求解，可得最大利润由3360元增加到3362元，其中增加的2元就是劳动时间的影子价格．因此，若可以聘用临时工人以增加劳动时间，付给临时工人的工资最多是每小时2元．其次，若还想知道以低于每小时2元的价格增加劳动时间，最多可购买多少劳动时间可以对目标函数以及关于劳动时间的约束条件作类似的修改，然后进行求解．例如，若以每小时元的价格聘用临时工人，最多可购买小时．⑶针对附加问题3，只需改变目标函数中的系数即可．将原来的目标函数改为：12max 9064z x x =+，约束条件不变．通过求解可得：最大利润有所增加，由原来的3360元增加到3720元，但生产计划没有改变，仍然是每天用20桶牛奶生产A ，用30桶牛奶生产B ． 4.1.4应用实例例1一个家庭有625英亩的土地可以用来种植农作物，这个家庭考虑种植的农作物有玉米、小麦和燕麦，预计可以有1000英亩-英尺的灌溉用水，农场工人每周可以投入的工作时间为300小时，其他的数据在表4-2中给出，为能够获得最大收益，每种作物应该种植多少表农场问题的有关数据玉 米 小 麦 燕 麦 现有量 灌溉用水（英亩-英尺） 劳力（人-小时/周） 收益（美元）4002002501000 300解 设应种植玉米1x 英亩，小麦2x 英亩和燕麦3x 英亩．可得如下线性规划模型：目标函数：123max 400200250z x x x =++约束条件：1231231233 1.510000.80.20.3300625x x x x x x x x x ++≤⎧⎪++≤⎨⎪++≤⎩1230,0,0x x x ≥≥≥使用Lingo 数学软件进行求解，编写程序如下：max =400*x1+200*x2+250*x3;x1+x2+x3<625; *x1+*x2+*x3<300; 3*x1+x2+*x3<1000; end程序运行结果为：应分别种植玉米187.5英亩，小麦437.5英亩和燕麦0英亩，获最大收益162500美元．例2 某市有甲、乙、丙、丁四个居民区，自来水由A 、B 、C 三个水库供应．四个区每天必须得到保证的基本生活用水量分别为30，70，10，10千吨，但由于水源紧张，三个水库每天最多只能供应50，60，50千吨自来水．由于地理位置的差别，自来水公式从各水库向各地区送水所需付出的引水管理费不同（见表4-3，其中C 水库与丁地区之间没有输水管道），其它管理费用都是450元／千吨．根据公司规定，各区用户按照统一标准900元／千吨收费．此外，四个区都向公司申请了额外用水量，分别为每天50，70，20，40千吨．该公司应如何分配供水量，才能获利最多为了增加供水量，自来水公司正在考虑进行水库改造，使三个水库每天的最大供水量都提高一倍，问那时供水方案应如何改变公司利润可增加多少表引水管理费引水管理费（元／千吨）甲 乙 丙 丁 A B C160 140 190130 130 200220 190 230170 150 ／解 决策变量应为A 、B 、C 三个水库(1,2,3)i =分别向甲、乙、丙、丁四个区(1,2,3,4)j =的供水量．设i 水库向j 区的日供水量为ij x 千吨，由于C 水库与丁区之间没有输水管道，即340x =，因此只有11个决策变量．于是可得如下线性规划模型： 目标函数：11121314max (450160)(450130)(450220)(450170)z x x x x =-+-+-+-21222324313233(450140)(450130)(450190)(450150)(450190)(450200)(450230)x x x x x x x +-+-+-+-+-+-+-约束条件：水库供应量限制可以表示为：1112131421222324313233506050x x x x x x x x x x x +++≤+++≤++≤ 各区基本用水量与额外用水量限制为：11213111213111213111213130807014010301050x x x x x x x x x x x x ≤++≤≤++≤≤++≤≤++≤0,1,2,3;1,2,3,4ij x i j ≥==使用Lingo 数学软件进行求解，编写程序如下：max =(450-160)*x11+(450-130)*x12+(450-220)*x13+(450-170)*x14+(450-140)*x21+(450-130)*x22+(450-190)*x23+(450-150)*x24 +(450-190)*x31+(450-200)*x32+(450-230)*x33; x11+x12+x13+x14<50; x21+x22+x23+x24<60; x31+x32+x33<50; x11+x21+x31<80; x11+x21+x31>30; x12+x22+x32<140; x12+x22+x32>70; x13+x23+x33<30; x13+x23+x33>10; x14+x24<50; x14+x24>10; end程序运行结果为：A 水库向乙区供水50千吨；B 水库向乙、丁区分别供水50，10千吨；C 水库向甲、丙分别供水40，10千吨．最大利润为47600元．对于本例来说，无论是目标函数还是约束条件都显得比较麻烦，特别是目标函数为多项相加．随着水库与居民区个数的增加，程序会更加复杂．下面来研究更一般的编程方法．为此，需假定C 水库向丁地区的引入管理费用为无穷大，本例可用1000元／千吨来代替．使用Lingo 数学软件中的高级编程技巧进行求解，编写程序如下：model:sets:sk/1..3/:g;dq/1..4/:sl1,sl2;link(sk,dq):c,x;endsetsdata:c=160 130 220 170140 130 190 150190 200 230 1000;g=50 60 50;sl1=30 70 10 10;sl2=80 140 30 50;enddata[obj] max=@sum(link(i,j):x(i,j)*(450-c(i,j)));@for(sk(i):@sum(dq(j):x(i,j))<g(i));@for(dq(j):@sum(sk(i):x(i,j))>sl1(j));@for(dq(j):@sum(sk(i):x(i,j))<sl2(j));end程序运行结果完全相同．如果三个水库每天的最大供水量都提高一倍，只需将数据段中的供水量修改成：g=100 120 100；或者对第一个约束条件作简单修改，在小于号后面将供水量扩大2倍，其它条件不变．最后的运行结果为：A水库向乙区供水100千吨；B水库向甲、乙、丁区分别供水30，40，50千吨；C水库向甲、丙分别供水50，30千吨．总利润为88700元．评注：本例考虑的是将某种物资从若干供应点运往一些需求点，在供需量约束条件下使总费用最少，或总利润最大．这类问题一般称为运输问题．注意：本例目标函数采用的是最大利润，而非最小成本．一般来说，成本最小，未必利润最大．当总收入是常数时，最小成本与最大利润是等价的；若总收入随决策变量的改变而变化时，最小成本与最大利润并不等价．通常追求的目标应该是最大利润，而非最小成本．非线性规划在工程技术、经济管理、交通运输和日常生活等诸多领域中，很多实际问题可以归结为线性规划问题，其目标函数和约束条件都是自变量（决策变量）的一次函数（线性函数）．但是，还有另外一些问题，其目标函数和（或）约束条件很难用线性函数表达．如果目标函数或约束条件中包含有非线性函数，就称这种规划问题为非线性规划问题．由于计算机技术的快速发展，使非线性规划的理论及其应用在近几十年来得以长驱进展．特别是Lingo数学软件的开发与应用，对非线性规划模型的求解提供了很大的帮助．4.2.1问题的提出1.引例1液体原料混合问题某公司将4种不同含硫量的液体原料（分别记为甲、乙、丙、丁）混合生产两种产品（分别记为A 、B ）．按照生产工艺的要求，原料甲、乙、丁必须首先倒入混合池中混合，混合后的液体再分别与原料丙混合生产A 和B ．已知原料甲、乙、丙、丁的含硫量分别是3，1，2，1（%），进货价格分别为6，16，10，15（千元／吨）；产品A 、B 的含硫量分别不能超过，（%），售价分别为9，15（千元／吨）．根据市场信息，原料甲、乙、丙的供应没有限制，原料丁的供应量最多为50吨；产品A 、B 的市场需求量分别为100、200（吨），问应如何安排生产2.引例2最佳选址问题某乡镇由12个主要的自然村组成，每个自然村的位置（用平面坐标x ，y 表示，距离单位：km ）和自然村的人口数（R ）如表4-4所示． 试根据需要解决如下问题：⑴ 目前准备在该乡镇建一个服务网点为各村提供各种服务，那么服务网点应该建在何处⑵ 假设各村人口增长了一倍，需要建两个服务网点，确定其位置．表最佳选址问题0 1234567 89101112X y R0 0 600 1000 800 1400 1200 700600800 1000 1200 1000 11004.2.2模型建立1.引例1液体原料混合问题的模型建立设11,y z 分别是产品A 中来自混合池和原料丙的吨数；22,y z 分别是产品B 中来自混合池和原料丙的吨数；混合池中原料甲、乙、丁所占的比例分别为124,,x x x ，优化目标是总利润（z ）最大．目标函数为：11221212412max 9()15()10()(61615)()z y z y z z z x x x y y =+++-+-+++1241124212(961615)(1561615)(910)(1510)x x x y x x x y z z =---+---+-+-约束条件为：⑴原料最大供应限制：412()50x y y +≤⑵产品最大需求限制：1122100,200y z y z +≤+≤ ⑶产品最大含硫量限制：12411124221122(3)2(3)22.5, 1.5x x x y z x x x y z y z y z ++++++≤≤++⑷其它限制：12412411221,,,,,,,0x x x x x x y z y z ++=≥2.引例2最佳选址问题的模型建立 （1）模型一的建立设服务网点的坐标为：(,)a b ；自然村的位置坐标：(,),1,2,,12i i x y i =；自然村的人口数：,1,2,,12i r i =，服务网点到各自然村的距离为：,1,2,,12i d i =．以自然村的人口数作为距离的权重，优化的目标为总距离最小．目标函数为：121min i ii z rd==∑约束条件为：1,2,,12i d i ==（2）模型二的建立设两个服务网点的坐标分别为：(,),1,2i i a b i =；自然村的位置坐标：(,)j j x y ，1,2,,12j =；自然村的人口数：,1,2,,12j r i =；服务网点i 到自然村j 的距离为：,1,2;1,2,,12ij d i j ==；服务网点i 对自然村j 服务的人口数为：,1,2ij c i =；1,2,,12j =；(),1,2k i i =表示第i 个服务网点服务的人口数占人口总数的比例．以服务网点对自然村服务的人口数作为距离的权重，优化的目标为总距离最小．目标函数为：12211min (,)(,)j i z c i j d i j ===⋅∑∑约束条件：12112121(,)(,)()()=()2()(,)2()j j i d i j c i j e i e i k i r j c i j r j ===⎧=⎪⎪=⎪⎪⎪⎨⎪⎪⎪⎪≥⎪⎩∑∑∑从以上两个引例可以总结出非线性规划的一般模型： 目标函数：12max(min)(,,,)n z f x x x =约束条件：1212(,,,)0,1,2,,(,,,)0,1,2,,i n j n h x x x i mg x x x j l ==⎧⎨≥=⎩目标函数为一般非线性函数，约束条件为一般非线性等式或非线性不等式．一般来说，目标函数与约束条件中只要有非线性项存在，即为非线性规划．特别地，若某非线性规划的目标函数为决策变量的二次函数，约束条件又全是线性的，就称这种规划为二次规划．4.2.3模型求解1.引例1液体原料混合问题的模型求解使用Lingo 数学软件进行求解，编写程序如下：max =(9-6*x1-16*x2-15*x4)*y1+(15-6*x1-16*x2-15*x4)*y2+(1-10)*z1+(15-10)*z2; x4*(y1+y2)<50;y1+z1<100;y2+z2<200;((3*x1+x2+x4)*y1+2*z1)/(y1+z1)<; ((3*x1+x2+x4)*y2+2*z2)/(y2+z2)<; x1+x2+x4=1; end用Lingo 求解时，如果怀疑不是全局最优解，用＂LINGO|OPTIONS ＂菜单命令启动＂Global Solver ＂选项卡上的＂Use Global Solver ＂选项，然后求解，可以得到全局最优解如下：24220.5,100x x y z ====，其余为0；目标函数值为450．如果将产品最大含硫量限制的约束条件作简单修改后，也可直接进行求解，并得到相同的结果．修改后的程序如下：max =(9-6*x1-16*x2-15*x4)*y1+(15-6*x1-16*x2-15*x4)*y2+(1-10)*z1+(15-10)*z2; x4*(y1+y2)<50;y1+z1<100;y2+z2<200;!((3*x1+x2+x4)*y1+2*z1)/(y1+z1)<; (3*x1+x2+**z1<0;!((3*x1+x2+x4)*y2+2*z2)/(y2+z2)<; (3*x1+x2+*y2+*z2<0; x1+x2+x4=1; end2.引例2最佳选址问题的模型求解针对模型一，使用Lingo 数学软件进行求解，编写程序如下： model :title :最佳选址（一）; sets :point/1..12/:x,y,r,dis; endsets data :X=0 ; Y=0 ;r=600 1000 800 1400 1200 700 600 800 1000 1200 1000 1100; enddata@for (point(i): dis(i)=((x(i)-a)^2+(y(i)-b)^2)^(1/2)); min =@sum (point: dis*r); end用Lingo 求解得到结果为： 3.601, 6.514a b ==．针对模型二，若取(1)(2)0.5k k==，使用Lingo数学软件进行求解，编写程序如下：model:title:最佳选址（一）;sets:point/1..12/:x,y,r;weizhi/1..2/:a,b,e;link(weizhi,point):c;endsetsdata:X=0 ;Y=0 ;r=600 1000 800 1400 1200 700 600 800 1000 1200 1000 1100;enddatasubmodel xuanzhi:min=@sum(link(i,j): c(i,j)*((a(i)-x(j))^2+(b(i)-y(j))^2)^(1/2));@for(weizhi(i): @sum(point(j):c(i,j))=e(i));@for(point(j): @sum(weizhi(i):c(i,j))) >2*r(j));endsubmodelcalc:e(1)=@sum(point:r);e(2)=@sum(point:r);@solve(xuanzhi);@ole('选址.xls','最佳位置a')=a;@ole('选址.xls','最佳位置b')=b;@ole('选址.xls','最优方案')=c;EndcalcEnd用Lingo求解得到结果为：两个服务网点的位置坐标为：(1.92,7.70)；(5.70,5.00)各服务网点服务人数对照表见表4-5．表服务人数对照表（限制服务网点的服务人数相同）1 2 3 4 5 6 7 8 9 10 11 12 人口总数（千人）网点1 网点2 02 00 020 2针对模型二，若不考虑服务网点服务人数的限制，使用Lingo数学软件进行求解，编写程序如下：model:title:最佳选址（一）;sets:point/1..12/:x,y,r;weizhi/1..2/:a,b,e;link(weizhi,point):c;endsets data :X=0 ; Y=0 ;r=600 1000 800 1400 1200 700 600 800 1000 1200 1000 1100; enddatasubmodel xuanzhi:min =@sum (link(i,j): c(i,j)*((a(i)-x(j))^2+(b(i)-y(j))^2)^(1/2)); !@for(weizhi(i): @sum(point(j):c(i,j))=e(i)); @for (point(j): @sum (weizhi(i):c(i,j))>2*r(j)); endsubmodel calc :!e(1)=@sum(point:r); !e(2)=@sum(point:r); @solve (xuanzhi);@ole ('选址.xls','最佳位置a')=a; @ole ('选址.xls','最佳位置b')=b; @ole ('选址.xls','最优方案')=c; endcalc用Lingo 求解得到结果为：两个服务网点的位置坐标为：(6.434,3.411)；(2.540,7.936)，各服务网点服务人数对照表见表4-6．表服务人数对照表（服务网点服务人数不限制）1 2 3 4 5 6 7 8 9 10 11 12 人口总数（千人）网点1 网点2 0 0 2 0 0 0 0 0 0 2 0 024.2.4应用实例例1求解二次规划：221212min 810z x x x x =+--12326x x +≤120,0x x ≥≥解 本例编写简单的Lingo 程序即可求解，编写Lingo 程序如下：max =8*x1+10*x2-x1^2-x2^2;3*x1+2*x2<6; 求解结果为：120.308, 2.538x x ==；目标函数值为：min 21.308z =例2 一个飞行管理问题在约10,000米高空的某边长160公里的正方形区域内，经常有若干架飞机作水平飞行．区域内每架飞机的位置和速度向量均由计算机记录其数据，以便进行飞行管理．当一架欲进入该区域的飞机到达区域边缘时，记录其数据后，要立即计算并判断是否会与区域内的飞机发生碰撞．如果会碰撞，则应计算如何调整各架（包括新进入的）飞机飞行的方向角，以避免碰撞．现假定条件如下：⑴不碰撞的标准为任意两架飞机的距离大于8公里；⑵飞机飞行方向角调整的幅度不应超过30度；⑶所有飞机飞行速度均为每小时800公里；⑷进入该区域的飞机在到达区域边缘时，与区域内飞机的距离应在60公里以上；⑸最多需考虑6架飞机；⑹不必考虑飞机离开此区域后的状况．请你对这个避免碰撞的飞行管理问题建立数学模型．列出计算步骤，对以下数据进行计算（方向角误差不超过度）．要求飞机飞行方向角调整的幅度尽量小．设该区域4个顶点的座标为(0,0)，(160,0)，(160,160)，(0,160)．记录数据为：飞机编号横座标x 纵座标y 方向角（度）1 150 140 2432 85 85 2363 150 1554 145 50 1595 130 150 230新进入 0 0 52注：方向角指飞行方向与x轴正向的夹角．试根据实际应用背景对你的模型进行评价与推广．解模型一：圆状模型下面是本例圆状模型建模的全过程．１.问题分析根据题目的条件，可将飞机飞行的空域视为二维xoy平面中的一个正方形区域，顶点(0,0)，(160,0)，(160,160)，(0,160)．各架飞机的飞行方向角为飞行方向与x 轴正向夹角（转角）．根据两飞机不碰撞的标准为二者距离大于8km，可将每架飞机视为一个以飞机坐标点为圆心、以4km为半径的圆状物体（每架飞机在空域中的状态由圆心的位置矢量和飞行速度矢量确定）．这样两架飞机是否碰撞就化为两圆在运行过程中是否相交的问题．两圆是否相交只要讨论它们的相对运动即可．2.模型假设⑴飞机进入区域边缘时，立即作出计算，每架飞机按照计算后的指示立即作方向角改变；⑵每架飞机在整个过程中至多改变一次方向；⑶忽略飞机转向的影响(转弯半径和转弯时间的影响)；⑷新飞机进入空域时，已在空域内部飞行的飞机的飞行方向已调合适，不会碰撞；⑸对每架飞机方向角的相同调整量的满意程度是一样的．3.模型的建立符号说明：i，j表示第i，第j架飞机的圆心；表示第i架飞机与第j架飞机的碰撞角，是两圆公切线的夹角中指向圆的那个ij角的一半，ij ji αα=；ij v 表示第i 架飞机相对于第j 架飞机的相对飞行速度； ij r 表示第i 架飞机与第j 架飞机的圆心距；ij β表示第i 架飞机对于第j 架飞机的相对速度与两架飞机圆心连线的夹角．规定以第i 架飞机为原点，i →j 连线从i 指向j 为正方向，逆时针旋转为正角，顺时针旋转为负角；AB ，CD 为两圆的公切线，//,//im AB in CD ． 另外再引入记号：i θ表示第i 架飞机的飞行方向与x 轴正向的夹角（转角）；(,)i i x y 表示第i 架飞机在坐标系中的位置矢量； i v 表示第i 架飞机的飞行速度矢量．由前面的分析将飞机作为圆状模型进行研究．两圆不相交，则表明不会发生碰撞事故；若两圆相交，则表明会发生碰撞事故．为了研究两飞机相撞问题，采用相对速度作为研究对象，因为飞机是否相撞的关键是相对速度，图4-1给出任意两架飞机间的关系．图4-1 第i 架飞机与第j 架飞机的碰撞角由图4-1中的关系得到两飞机不相撞（两圆不相交）的充要条件是||||ij ij βα>．当||||ij ij βα≤时，则通过调整两架飞机的方向角i θ，j θ，使飞机不相撞．首先讨论相对飞行速度的方向角ij β的改变量ij β∆与第i 、第j 架飞机方向角的改变量i θ∆， j θ∆的关系．由题目条件知，对于第i 架飞机：||800i v km A ==．设第i ，j 架飞机改变飞行方向前的速度分别为：i 1i i v Ae θ=，i 1j j v Ae θ=；改变飞行方向后的速度分别为：i()2i i i v Ae θθ+∆=，i()2j j j v Aeθθ+∆= ．则飞行方向改变前后的相对速度分别为：i i 111()j i ij i j v v v A e e θθ=-=- i()i()222()j j i i ij i j v v v A e e θθθθ+∆+∆=-=-2i()i()i 1i ()()j j i i j i ijij v A e e v A e e θθθθθθ+∆+∆-=-cos()isin()cos()isin()cos isin cos isin i i i i j j j j i j j jθθθθθθθθθθθθ+∆++∆-+∆-+∆=+--2sin(sin i cos)222sin(sini cos)222i i j ji i j ji i j ji j i ji jθθθθθθθθθθθθθθθθθθ+∆--∆+∆++∆+∆++∆-=-++-i2sin 2sin 2i ji i j ji jeθθθθθθθθ∆+∆+∆--∆=-即2ij v 与1ijv 辐角相差2i jθθ∆+∆．将其归纳为：定理 对第i ，j 架飞机，其相对速度方向角的改变量ij β∆等于两飞机飞行方向角的改变量之和的一半2i jθθ∆+∆．由题目的要求调整飞行方向角时不能超过30°，即|i θ∆|≤30 , i=1,2,…,6． 要保证调整飞行方向后飞机不碰撞，应有： ||ij ij ij ββα+∆≥ 于是可得如下非线性规划模型： 目标函数：621min ||i i θ=∆∑约束条件：||,,1,2,,6,2||301,2,,6ij ij ij i j ij i i j i j i ββαθθβθ+∆≥=≠⎧⎪∆+∆⎪∆=⎨⎪∆≤=⎪⎩， 其中ij β，ij α可由题中已知的参数计算得到：=arcsin(8/)ij ji ij r αα=；ij r =22sin sin arctan arctan cos cos i j j i ij i jj iy y x x θθβθθ--=--（）-（）； 其中，2arctan b a （）与arctan b a （）的区别为：2arctan ba（）表示取值位于π-到π之间的辐角：可根据点(,)a b 所在的象限确定．由此计算得到的ij β取值位于2π-到2π之间，还需要将它转换到π-到π之间（超过π时就减去2π；小于π-就加上2π）．4.模型求解计算任两架飞机间的参数ij β，编写Matlab 指令如下： clear,clcx=[ 150 140 243;85 85 236;150 155 ;... 145 50 159;130 150 230;0 0 52]; x0=x(:,1); y0=x(:,2);xita=x(:,3)*pi/180; b=zeros(6); for i=1:6for j=i+1:6x11=cos(xita(i))-cos(xita(j)); x12=sin(xita(i))-sin(xita(j)); if x11>=0b(i,j)=atan(x12/x11); elseif x12>=0b(i,j)=pi+atan(x12/x11); elseif x12<0b(i,j)=-pi+atan(x12/x11); endx21=x0(j)-x0(i); x22=y0(j)-y0(i); if x21>=0b(i,j)=b(i,j)-atan(x22/x21); elseif x22>=0b(i,j)=b(i,j)-atan(x22/x21)-pi;elseif x22<0b(i,j)=b(i,j)-atan(x22/x21)+pi;endif b(i,j)>pib(i,j)=b(i,j)-2*pi;elseif b(i,j)<-pib(i,j)=b(i,j)+2*pi;endendendb=b*180/pi;save beta b计算结果见表4-7．β的值表使用Matlab计算得到的ij1 2 3 4 5 612３４５６α可以在Lingo中直接计算．编写求解该问题的Lingo程序如下：对于ijmodel:title:飞行管理问题的非线性规划模型一;sets:plane/1..6/:x0,y0,d_cita;! d_cita为调整的角度;link(plane,plane)|&1#lt#&2:alpha,beta;endsetsdata:x0,y0=150 14085 85150 155145 50130 1500 0 ;beta=; enddata !计算alpha;@for (link(i,j):@sin (alpha*3./180)=8/((x0(i)-x0(j))^2+(y0(i)-y0(j))^2)^.5); @for (link(i,j):@abs (beta(i,j)+*d_cita(i)+*d_cita(j))>alpha(i,j);); @for (link:@bnd (0,alpha,90));@for (plane:@bnd (-30,d_cita,30)); [obj]min =@sum (plane:(d_cita)^2); end5.结果检验对题目所给实例进行计算得如下调整方案：10θ∆=， 20θ∆=， 3 2.062465θ∆=， 4-0.4954514θ∆=， 50θ∆=， 6 1.567013θ∆=．各飞行方向角按此方案调整后，系统各架飞机均满足||ij ij βα>（即不会相撞）．其中有些飞机对可能会有||0.01ij ij βα-<（°是题目要求的计算精度）．如果希望||0.01ij ij βα≥+，只须将模型中的ij α用0.01ij ij αα=+代替即可．6.模型评价与改进此模型采用圆状模型分析碰撞问题是合理的，同时采用相对速度作为判断标准，既体现了碰撞的本质（相对运动），又简化了模型的计算．题目要求飞机飞行方向角调整的幅度尽量小，这个尽量小是针对每架飞机而言的，同时也要求整体满意程度（即对管理层而言，应使每架飞机的调整都尽量的小）．因此构造目标函数时，也可以认为若对方向角调整量最大的飞机而言，其调整量可满意，则由假设（5）对其余飞机调整量均可满意．即要求每架飞机的调整量都小于某个数θ(0)θ≥．故目标函数也可取：min θ．于是可得如下线性规划模型：目标函数：min θ约束条件：||,,1,2,,6,2||301,2,,6||1,2,,6ij ij ij i j iji i i j i j i i ββαθθβθθθ+∆≥=≠⎧⎪∆+∆⎪∆=⎪⎨⎪∆≤=⎪∆≤=⎪⎩，， 模型二： 最短距离模型１.问题分析目标函数的选取与模型一相同．进入该区域的飞机在到达该区域边缘时，与区域内的飞机的距离应在60km 以内，很容易验证目前所给数据是满足的，因此，约束条件只需要限制任意两架位于该区域内的飞机的距离应大于8km ．但这个问题的难点在于飞机是动态的，这个约束不好直接描述．为此，可以考虑在飞行过程中任意两架飞机的最短距离大于8km 即可．飞行时间可以只考虑一架飞机飞越该区域所需的最长时间，若超过这个时间，即使两架飞机的最短距离小于8km ，由于飞机已经离开该区域，因此不再予以考虑．2.模型假设 与模型一相同．3.模型的建立 符号说明：i θ表示第i 架飞机的飞行方向与x 轴正向的夹角（转角）；00(,)i i x y 表示第i 架飞机在调整前的位置坐标； (,)t t i i x y 表示第i 架飞机t 时刻的位置坐标； t 表示飞机的飞行时间； v 表示飞机的飞行速度；i T 表示第i 架飞机飞出区域的时刻；max T 表示任意一架飞机在该区域内停留的最长时间； min{,}ij i j T T T =；*ijt 表示第i 架飞机与第j 架飞机距离最短的时刻； i θ表示第i 架飞机的飞行方ij r 表示第i 架飞机与第j 架飞机的距离；2()[]64ij ij f t r =-记飞机飞行速率为v ，以当前时刻为0时刻，设第i 架飞机在调整前的位置坐标为00(,)i i x y ，t 时刻的位置坐标(,)t t i i x y ，即00cos ,sin t t i i i i i i x x vt y y vt θθ=+=+如果要严格表示两架位于该区域内的飞机的距离应大于8km ，则需要考虑每架飞机在区域内到为飞行时间的长度．记i T 为第i 架飞机飞出区域的时刻，即00min{0:cos 0,160;sin 0,160}i i i i i T t x vt y vt θθ=>+=+=记t 时刻第i 架飞机与第j 架飞机的距离为ij r ，并记2()[]64ij ij f t r =-，这时在区域内飞机不相撞的约束条件就变成了2()[]640,(0)ij ij ij f t r t T =-≥≤≤ 其中min{,}ij i j T T T = 此外，经过计算可以得到：002002()(cos cos )(sin sin )64ij i i j j i i j j f t x vt x vt y vt y vt θθθθ=+--++---2ij ij ij ij z b z c =++其中2sin2i jij z vt θθ-=00002[()sin()cos22i ji jij i j i j b x x y y θθθθ++=--+-002002()()64ij i j i j c x x y y =-+--所以，()ij f t 是一个关于ij z 的二次函数，当2ij ij b z =-，即/4sin2i jij t b v θθ-=-（记为*ij t ）时，函数()ij f t 取最小值2/4ij ij b c -+．若*0ij t <，只要初始时刻不相撞即可，此时满足条件，不需要限制； 若*ij ij t T ≥，只需要()0ij ij f T ≥即可；若*0ij ij t T <<，*2()/40ij ij ij ij f t b c =-+≥即可，即．实际上，()0ij ij f T ≥在*0ij ij t T <<时也成立，因此，可以不再附加*ij ij t T ≥的条件，于是可得如下非线性规划模型：目标函数：621min ||i i θ=∆∑约束条件：2*||301,2,,6()040,(0)i ij ij ijij ij ij i f T b c t T θ⎧∆≤=⎪≥⎨⎪-≤<<⎩，4.模型求解由于ij T 的计算相当复杂，求解时可进一步简化：不单独考虑每架飞机在区域内停留的时间，而以最大时间max 0.283()T h ==代替，此时所有max =ij T T ．实际上强化了问题的要求，即考虑了有些飞机可能已经飞出区域，但仍不允许两架飞机的距离小于8km ．这个简化的模型可以如下输入Lingo软件：model:title: 飞行管理问题的非线性规划模型二;sets:plane/1..6/:x0,y0,cita0,cita1,d_cita;!cita0表示初始角度，cita1为调整后的角度，d_cita为调整的角度;link(plane,plane)|&1#lt#&2:b,c;endsetsdata:x0,y0,cita0=150 140 24385 85 236150 155145 50 159130 150 2300 0 52;max_cita=30;t_max=;v=800;pi=3.;enddatainit:d_cita=0 0 0 0 0 0;endinit@for(plane:cita1-cita0=d_cita);@for(link(i,j):b(i,j)=-2*(x0(i)-x0(j))*@sin((cita1(i)+cita1(j))*pi/360)+2*(y0(i)-y0(j))*@cos((cita1(i)+cita1(j))*pi/360);c(i,j)=(x0(i)-x0(j))^2+(y0(i)-y0(j))^2-64;);!避免碰撞的条件;!右端点非负;@for(link(i,j):[right](2*v*t_max*@sin((cita1(i)-cita1(j))*pi/360))^2+b(i,j)*(2*v*t_max*@sin((cita1(i)-cita1(j))*pi/360))+c(i,j)>0);!左端点非负;@for(link(i,j):c(i,j)>0);!最小点非负;@for(link(i,j):[minimum]@if(-b(i,j)/4/v/@sin((cita1(i)-cita1(j))*pi/360)#g t#0#and#-b(i,j)/4/v/@sin((cita1(i)-cita1(j))*pi/360)#lt#t_max,b(i,j)^2-4*c(i,j),-1)<0);!@for(link(i,j):b(i,j)^2-4*c(i,j)<0);@for(link:@free(b));。