数学建模竞赛案例选讲

数学建模竞赛案例选讲,飞行管理问题,1995年A题,飞行管理问题（1995年全国大学生数学建模竞赛试题A）,问题：在约10000米高空的正方形区域内，有若干架飞机作水平飞行。区域内每架飞机的位置和速度向量均由计算机记录数据，以便进行飞行管理，当一架欲进入该区域的飞机到达区域边缘时，记录其数据后，要立即计算并判别是否会与区域内的飞机发生碰撞。如果会碰撞，则应计算如何调整各架（包括新进入的）飞机飞行的方向角，以避免碰撞。,假定条件:,1.不碰撞的标准是任意两架飞机的距离大于8km；,2.飞机飞行方向角调整幅度不应超过30度，而要尽可能小；,3.所有飞机的飞行速度为800km/h，不受其他因素影响；,4.进入该区域的飞机在到达边缘时，与该区域内的飞机的距离应在60km以上；,5.不考虑飞机离开区域后的情况；,6.建模时暂考虑6架飞机；,问题的提出：,请你对这个避免碰撞的飞行管理问题建立数学模型，列出计算步骤，对以下数据进行计算（方向角误差不超过0.01度），要求飞机的方向角调整的幅度尽量小。,设该区域4个顶点的坐标为(0，0),(160，0),(160，160),(0，160），记录数据如下表（其中方向角指飞行方向与x轴正向夹角),问题分析,根据题目的条件，可将飞机飞行的空域视为二维平面xoy中的一个正方形，顶点在(0，0),(160，0),(160，160),(0，160）。各架飞机的飞行方向角为飞行方向与x轴正向夹角（转角）。根据两飞机不碰撞的标准为二者距离大于8km，可将每架飞机视为一个以飞机为圆心、以4为半径的圆状物体（每架飞机在空域中的状态由圆心的位置矢量和飞行速度矢量确定）。这样两架飞机是否碰撞就化为两圆在运动中是否相交的问题。两圆是否相交只要讨论它们的相对运动即可。,C,建模时补充假定条件:,1.飞机在所定区域内作直线飞行，不偏离航向；,2.飞机管理系统内不发生意外，如发动机失灵，或其他意外原因迫使飞机改变航向；,4.飞机管理系统发出的指令应被飞机立即执行，即认为转向是瞬间完成的（忽略飞机转向的影响，即转弯半径和转弯时间的影响）；,3.飞机进入区域边缘时，立即作出计算，每架飞机按照计算后的指示立即作方向角改变；,5.每架飞机在在整个过程中指点改变一次方向,6.新飞机进入区域时，已在区域内部的飞机的飞行方向已调整合适，不会碰撞；,7.对每架飞机方向角的相同调整量的满意程度是一样的。,模型的建立,（1）圆状模型,采用相对速度作为研究对象,符号说明：,i,j第i，第j架飞机的圆心；,ij第i，第j架飞机的碰撞角，ij=ji;,vij第i架飞机相对第j架飞机的相对飞行速度;,lij第i，第j架飞机的圆心距,i第i架飞机的飞行方向与x轴正向夹角（逆时针为正）,xi第i架飞机的位置矢量,vi第i架飞机的的速度矢量,ij第i飞机对第j架飞机的相对速度与两架飞机圆心连线的夹角（逆时针为正）,不碰撞,|ij|ij,(2)由圆状模型导出的方程,讨论ij的改变量与第i第j两架飞机飞行方向角改变量i,j的关系,由题目条件知|vi|=A=800，可用复数表示速度,设第i，j飞机飞行方向改变前的速度分别为,改变后的速度分别为,改变前后相对速度分别为,两者之商的幅角就是ij,定理：对第i，第j两架飞机，其相对速度方向ij的改变量ij等于两飞机飞行方向角改变量之和的一半，即,模型,目标函数：,Min其中为各飞机方向角调整量的最大值,或为,约束条件：,调整方向角时不能超过300：,调整飞行方向后飞机不能碰撞：,模型为,或为,化为线性规划模型,由于i可正可负，为使各变量均非负，引入新变量：,模型化为,模型求解,ij的计算,model:sets:plane/1.6/:x0,y0;link(plane,plane):alpha,sin2;endsetsfor(link(i,j)|i#ne#j:sin2(i,j)=64/(x0(i)-x0(j)2+(y0(i)-y0(j)2););for(link(i,j)|i#ne#j:(sin(alpha*3.14159265/180.0)2=sin2;);data:x0=150,85,150,145,130,0;y0=140,85,155,50,150,0;enddataend,ALPHA(1,1)1.234568ALPHA(1,2)5.391190ALPHA(1,3)752.2310ALPHA(1,4)5.091816ALPHA(1,5)2000.963ALPHA(1,6)2.234507ALPHA(2,1)5.391190ALPHA(2,2)1.234568ALPHA(2,3)4.804024ALPHA(2,4)6.613460ALPHA(2,5)5.807866ALPHA(2,6)3.815925ALPHA(3,1)752.2310ALPHA(3,2)4.804024ALPHA(3,3)1.234568ALPHA(3,4)4.364672ALPHA(3,5)1102.834ALPHA(3,6)2.125539ALPHA(4,1)5.091816,ALPHA(4,2)6.613460ALPHA(4,3)4.364672ALPHA(4,4)1.234568ALPHA(4,5)4.537692ALPHA(4,6)2.989819ALPHA(5,1)2000.963ALPHA(5,2)5.807866ALPHA(5,3)1102.834ALPHA(5,4)4.537692ALPHA(5,5)1.234568ALPHA(5,6)2.309841ALPHA(6,1)2.234507ALPHA(6,2)3.815925ALPHA(6,3)2.125539ALPHA(6,4)2.989819ALPHA(6,5)2.309841ALPHA(6,6)1.234568,整理可得ij的值（单位角度）,也可以用MATLAB计算ij的值,x=150,85,150,145,130,0;y=140,85,155,50,150,0;k=length(x);alpha=zeros(k);fori=1:kforj=1:kifi=jalpha(i,j)=0;elsealpha(i,j)=(180/3.14159265)*asin(8/sqrt(x(i)-x(j)2+(y(i)-y(j)2);endendendalpha,计算ij的值程序为,计算结果为,05.3911902372235.391190237223032.2309526723314.8040239337975.0918164485506.61346048987220.9633608931285.8078662434212.2345067369953.815924775399,32.2309526723315.0918164485504.8040239337986.61346048987204.3646718991114.364671899111022.8336542040094.5376924624022.1255388575512.989819139045,20.9633608931282.2345067369955.8078662434213.81592477539922.8336542040092.1255388575514.5376924624032.98981913904502.3098413654052.3098413654050,ij的计算:,a=150,85,150,145,130,0;b=140,85,155,50,150,0;x=a+b*i;c=243,236,220.5,159,230,52*pi/180;v=exp(i*c);k=length(a);fori=1:kforj=1:kbeita(i,j)=(angle(v(i)-v(j)-angle(x(j)-x(i)*180/pi;endendbeita,用matlab程序编写,beita=0109.2636-128.250024.1798-186.934914.4749109.26360-88.8711-42.2436-92.30489.0000231.7500271.1289012.4763301.21380.310824.1798-42.243612.476305.9692-3.5256173.0651267.6952-58.78625.969201.914414.47499.00000.3108-3.52561.91440,运算结果,最优解的计算,用LINGO求解,程序如下,model:sets:plane/1.6/:cita;link(plane,plane):alpha,beta;endsetsmin=sum(plane:abs(cita);for(plane(i):bnd(-30,cita(i),30););for(link(i,j)|i#ne#j:abs(beta(i,j)+0.5*cita(i)+0.5*cita(j)alpha(i,j););,data:alpha=0.000000,5.391190,32.230953,5.091816,20.963361,2.234507,5.391190,0.000000,4.8040024,6.813460,5.807866,3.815925,32.230953,4.804024,0.000000,4.364672,22.833654,2.125539,5.091816,6.613460,4.363673,0.000000,4.537692,2.989819,20.963361,5.807866,22.833654,4.537692,0.000000,2.309841,2.234507,3.815925,2.125539,2.989819,2.309841,0.000000;beta=0.000000109.263642-128.25000024.179830173.06505113.474934109.2636420.000000-88.871096-42.243563-92.3048479.000000-128.250000-88.870960.00000012.476311-58.7862430.31080924.179830-42.24356312.4763110.0000005.969234-3.525606174.065051-92.304846-58.7862445.9692340.0000001.91438314.4749349.0000000.310809-3.5256061.9133830.000000;enddataend,用MATLAB计算编程如下,functionf,g=plane(x)alpha=0.000000,5.391190,32.230953,5.091816,20.963361,2.234507,5.391190,0.000000,4.8040024,6.813460,5.807866,3.815925,32.230953,4.804024,0.000000,4.364672,22.833654,2.125539,5.091816,6.613460,4.363673,0.000000,4.537692,2.989819,20.963361,5.807866,22.833654,4.537692,0.000000,2.309841,2.234507,3.815925,2.125539,2.989819,2.309841,0.000000;beta=0.000000109.263642-128.25000024.179830173.06505113.474934109.2636420.000000-88.871096-42.243563-92.3048479.000000-128.250000-88.870960.00000012.476311-58.7862430.31080924.179830-42.24356312.4763110.0000005.969234-3.525606174.065051-92.304846-58.7862445.9692340.0000001.91438314.4749349.0000000.310809-3.5256061.9133830.000000;f=abs(x(1)+abs(x(2)+abs(x(3)+abs(x(4)+abs(x(5)+abs(x(6);g(1)=alpha(1,2)-abs(beta(1,2)+0.5*x(1)+0.5*x(2);g(2)=alpha(1,3)-abs(beta(1,3)+0.5*x(1)+0.5*x(3);g(3)=alpha(1,4)-abs(beta(1,4)+0.5*x(1)+0.5*x(4);g(4)=alpha(1,5)-abs(beta(1,5)+0.5*x(1)+0.5*x(5);,g(5)=alpha(1,6)-abs(beta(1,6)+0.5*x(1)+0.5*x(6);g(6)=alpha(2,3)-abs(beta(2,3)+0.5*x(2)+0.5*x(3);g(7)=alpha(2,4)-abs(beta(2,4)+0.5*x(2)+0.5*x(4);g(8)=alpha(2,5)-abs(beta(2,5)+0.5*x(2)+0.5*x(5);g(9)=alpha(2,6)-abs(beta(2,6)+0.5*x(2)+0.5*x(6);g(10)=alpha(3,4)-abs(beta(3,4)+0.5*x(3)+0.5*x(4);g(11)=alpha(3,5)-abs(beta(3,5)+0.5*x(3)+0.5*x(5);g(12)=alpha(3,6)-abs(beta(3,6)+0.5*x(3)+0.5*x(6);g(13)=alpha(4,5)-abs(beta(4,5)+0.5*x(4)+0.5*x(5);g(14)=alpha(4,6)-abs(beta(4,6)+0.5*x(4)+0.5*x(6);g(15)=alpha(5,6)-abs(beta(5,6)+0.5*x(5)+0.5*x(6);,执行程序,x0=0,0,0,0,0,0;v1=-30*ones(1,6);v2=30*ones(1,6);opt=;x=constr(plane,x0,opt,v1,v2),结果：,x=-0.00000576637983-0.000005766379832.58794980234726-0.000012434879850.000036204730951.04151019765274,最优解：,模型检验,各飞行方向按此方案调整后，系统各架飞机均满足|ij+(I+j)/2|ij,结果是正确的。,模型的评价与推广,（1）此模型采用圆状模型分析碰撞问题是合理的，同时采用相对速度作为判别标准，既体现了碰撞的本质（相对运动），又简化了模型的计算；,（2）建模中用了适当的简化，将一个复杂的非线性规划问题简化为线性规划问题，既求到合理的解，又提高了运算速度，这对解决高速飞行的飞机碰撞问题是十分重要的。此模型对题目所提供的例子计算得出的结果是令人满意的。,（3）由对称性知模型中的约束个数是（n是飞机数），所有约束条件数是，计算量增加不大。,投资的收益和风险,1998年A题,问题的提出：,市场上有n种资产（如股票、债券、）Si(i=1,n)供投资者选择，某公司有数额为M的一笔相当大的资金可用作一个时期的投资。公司财务分析人员对这n种资产进行了评估，估算出在这一时期内购买Si的平均收益率为，并预测出购买Si的风险损失率为。考虑到投资越分散，总的风险越小，公司确定，当用这笔资金购买若干种资产时，总体风险可用所投资的Si中最大的一个风险来度量。,购买Si要付交易费，费率为，并且当购买额不超过给定值时，交易费按购买计算（不买当然无须付费）。另外，假定同期银行存款利率是,且既无交易费又无风险。（=5%）,已知n=4时的相关数据如下：,试给该公司设计一种投资组合方案，即用给定的资金，有选择地购买若干种资产或存银行生息，使净收益尽可能大，而总体风险尽可能小。,2.试就一般情况对以上问题进行讨论，并利用以下数据进行计算。,假定条件:,1.题中所给利率均为一年，投资期为一年。2.公司的资金足够多，且全部用于投资或存入银行。3.投资期内不再做其他交易，利润仅在期末实现。4.风险损失率指投资到期后，如果风险发生，损失额占投资额的百分比。5.对于风险出现的概率与收益的波动，在本模型中不予考虑。6.总体的风险用最大的风险来衡量。7.当购买额不超过给定值ui时，交易费按购买ui计算。,符号说明:,M：投资总额Xi：对第i种投资项目的投资（不含交易费）占总投资额的比例,Q：总体风险ri：对第i种投资项目的平均收益率qi：对第i种投资项目的风险损失率pi：对第i种投资项目的交易费率ui：对第i种投资项目的交易费应付最小计算额ei：对第i种投资项目的交易费：乐观度常数,模型的建立:,考虑两个主要问题：,1.投资所获得的总收益尽量大2.投资所承担的风险尽量小,不考虑收益波动,用xiri度量收益,不考虑出现风险的概率,以当风险发生时的损失xiqi度量风险,数学模型:,这是一个多目标非线性规划,求解困难，需要简化,化为单目标规划,化为线性规划,化为单目标规划：,若有m个目标fi(x)，分别给以权系数i(i=1,2,m),然后作新目标函数（也称效用函数）,本问题,将非线性规划转化为线性规划规划：,目标函数线性化：,令,加入约束,约束条件线性化：,当投资额相当大时，可以认为xiui，,由于M是常数，可以从目标函数中去除,线性规划模型,即,模型求解,用MATLAB求解,计算线性规划时用命令:y=lp(C,A,b,v1,v2,x0,n),functionx=tzyh(S,a)m=length(a);n=length(S(:,1);h=1;v1=zeros(n+1,1);x0=eye(1,n+1);b=eye(n,1);fori=1:mforj=1:nC(j)=-a(i)*(1+S(j,1);endC(n+1)=1-a(i);D=ones(1,n)+S(:,3);D=D,0;K=S(2:n,2:2);,F=diag(K);H=zeros(n-1,1),F,-ones(n-1,1);A=D;H;y=lp(C,A,b,v1,x0,h);x(1,i)=a(i);forj=1:n+1x(1+j,i)=y(j);endendx,z=S(:,1:1)*x(2:n+1,:)fori=1:mforj=1:nQ(j,i)=x(j+1,i)*S(j,2);endendQ；f=x(n+2:n+2,:)；plot(f,z,b*)；,建立函数M文件tzyh如下：,其中S是由收益率，风险损失率，交易费率构成的矩阵，a是乐观程度参数数组。,S=0.05,0.28,0.21,0.23,0.25;0,0.025,0.015,0.055,0.0260,0.01,0.02,0.045,0.065S=S,问题1）的运算结果：,S=0.0500000.28000.02500.01000.21000.01500.02000.23000.05500.04500.25000.02600.0650,a=0:0.1:1,a=00.10000.20000.30000.40000.50000.60000.70000.80000.90001.0000,运行tzyh(S,a)得,最优解、收益、及最大风险情况,在最优解处各项投资的风险情况,对于固定的a，最优投资方案中各项投资的风险相同,又如：对a=0.077，,x=0.00000.23760.39600.10800.2284,收益：0.2316,风险：0.0059,问题2）的求解：,a=0:0.1:1,a=00.10000.20000.30000.40000.50000.60000.70000.80000.90001.0000,S=0.05,0.096,0.185,0.494,0.239,0.081,0.14,0.407,0.312,0.336,0.368,0.118,0.09,0.35,0.094,0.150,0.42,0.54,0.60,0.42,0.012,0.39,0.68,0.334,0.533,0.40,0.31,0.055,0.46,0.053,0.230,0.021,0.032,0.060,0.015,0.076,0.034,0.056,0.031,0.027,0.029,0.051,0.057,0.027,0.045,0.076S=S,运行tzyh(S,a)得,运行tzyh(S,a)得,S=0.0500000.09600.42000.02100.18500.54000.03200.49400.60000.06000.23900.42000.01500.08100.01200.07600.14000.39000.03400.40700.68000.05600.31200.33400.03100.33600.53300.0270,0.36800.40000.02900.11800.31000.05100.09000.05500.05700.35000.46000.02700.09400.05300.04500.15000.23000.0760,各项投资的风险情况：,数据分析,通过上述数据，我们发现如下特点：,（1）定理：,如果一个投资方案是所谓的最优非劣解，则各项投资,所冒的风险应基本