数学建模竞赛案例选讲.ppt

数学建模竞赛案例选讲,飞行管理问题,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.234568 ALPHA(1,2)5.391190 ALPHA(1,3)752.2310 ALPHA(1,4)5.091816 ALPHA(1,5)2000.963 ALPHA(1,6)2.234507 ALPHA(2,1)5.391190 ALPHA(2,2)1.234568 ALPHA(2,3)4.804024 ALPHA(2,4)6.613460 ALPHA(2,5)5.807866 ALPHA(2,6)3.815925 ALPHA(3,1)752.2310 ALPHA(3,2)4.804024 ALPHA(3,3)1.234568 ALPHA(3,4)4.364672 ALPHA(3,5)1102.834 ALPHA(3,6)2.125539 ALPHA(4,1)5.091816,ALPHA(4,2)6.613460 ALPHA(4,3)4.364672 ALPHA(4,4)1.234568 ALPHA(4,5)4.537692 ALPHA(4,6)2.989819 ALPHA(5,1)2000.963 ALPHA(5,2)5.807866 ALPHA(5,3)1102.834 ALPHA(5,4)4.537692 ALPHA(5,5)1.234568 ALPHA(5,6)2.309841 ALPHA(6,1)2.234507 ALPHA(6,2)3.815925 ALPHA(6,3)2.125539 ALPHA(6,4)2.989819 ALPHA(6,5)2.309841 ALPHA(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);for i=1:k for j=1:k if i=j alpha(i,j)=0;else alpha(i,j)=(180/3.14159265)*asin(8/sqrt(x(i)-x(j)2+(y(i)-y(j)2);end endendalpha,计算ij的值程序为,计算结果为,0 5.391190237223 5.391190237223 032.230952672331 4.80402393379720.963360893128 5.807866243421,32.230952672331 5.091816448550 0 4.364671899111 4.364671899111 0 22.833654204009 4.537692462402,20.963360893128 2.234506736995 4.537692462403 2.989819139045 0 2.309841365405 2.309841365405 0,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);for i=1:k for j=1:k beita(i,j)=(angle(v(i)-v(j)-angle(x(j)-x(i)*180/pi;endendbeita,用matlab程序编写,beita=0 109.2636-128.2500 24.1798-186.9349 14.4749 109.2636 0-88.8711-42.2436-92.3048 9.0000 231.7500 271.1289 0 12.4763 301.2138 0.3108 24.1798-42.2436 12.4763 0 5.9692-3.5256 173.0651 267.6952-58.7862 5.9692 0 1.9144 14.4749 9.0000 0.3108-3.5256 1.9144 0,运算结果,最优解的计算,用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.000000 109.263642-128.250000 24.179830 173.065051 13.474934109.263642 0.000000-88.871096-42.243563-92.304847 9.000000-128.250000-88.87096 0.000000 12.476311-58.786243 0.31080924.179830-42.243563 12.476311 0.000000 5.969234-3.525606174.065051-92.304846-58.786244 5.969234 0.000000 1.91438314.474934 9.000000 0.310809-3.525606 1.913383 0.000000;enddataend,用MATLAB计算编程如下,function f,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.000000 109.263642-128.250000 24.179830 173.065051 13.474934109.263642 0.000000-88.871096-42.243563-92.304847 9.000000-128.250000-88.87096 0.000000 12.476311-58.786243 0.31080924.179830-42.243563 12.476311 0.000000 5.969234-3.525606174.065051-92.304846-58.786244 5.969234 0.000000 1.91438314.474934 9.000000 0.310809-3.525606 1.913383 0.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.00000576637983 2.58794980234726-0.00001243487985 0.00003620473095,最优解：,模型检验,各飞行方向按此方案调整后，系统各架飞机均满足|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),function x=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);for i=1:m for j=1:n C(j)=-a(i)*(1+S(j,1);end C(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);for j=1:n+1 x(1+j,i)=y(j);endendx,z=S(:,1:1)*x(2:n+1,:)for i=1:mfor j=1:n Q(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.0500 0 0 0.2800 0.0250 0.0100 0.2100 0.0150 0.0200 0.2300 0.0550 0.0450 0.2500 0.0260 0.0650,a=0:0.1:1,a=0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000,运行tzyh(S,a)得,最优解、收益、及最大风险情况,在最优解处各项投资的风险情况,对于固定的a，最优投资方案中各项投资的风险相同,又如：对a=0.077，,x=0.0000 0.2376 0.3960 0.1080 0.2284,收益：0.2316,风险：0.0059,问题2）的求解：,a=0:0.1:1,a=0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.70000.8000 0.9000 1.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.15 0,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.23 0,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.0500 0 0 0.0960 0.4200 0.0210 0.1850 0.5400 0.0320 0.4940 0.6000 0.0600 0.2390 0.4200 0.0150 0.0810 0.0120 0.0760 0.1400 0.3900 0.0340 0.4070 0.6800 0.0560 0.3120 0.3340 0.0310 0.3360 0.5330 0.0270,0.3680 0.4000 0.0290 0.1180 0.3100 0.0510 0.0900 0.0550 0.0570 0.3500 0.4600 0.0270 0.0940 0.0530 0.0450 0.1500 0.2300 0.0760,各项投资的风险情况：,数据分析,通过上述数据，我们发现如下特点：,（1）定理：,如果一个投资方案是所谓的最优非劣解，则各项投资,所冒的风险应基本相等，即,引理：,在大资金的前提下，ui的问题可近似不考虑。,定理的证明：,用反证法。假设在最优方案x1,x2,xn中,有x1q1x2q2,为了简化问题，我们只考虑n=2，即两种方案时的证明,必小于或等于最优解处函数目标函数值,即,要是y正、负均成立，必须y系数为0，即,但此系数一般不为0，故此不等式难以成立。因此假设非真。,同理可证x1q1x2q2情况。故x1q1x2q2。可将两种方案推广到多种方案情况。定理得证。,（2）收益于风险的关系,在的变化中，随着的增加，收益在不断增加，风险也在增加。说明谨慎程度越强，风险越小但受益也越小。具有明显的实际意义，投资者可根据其对风险的承受能力和投资环境，选择适当的，然后在求出最优投资方案。,对于问题（1），我们用在01内的200个等分点，求出最优投资中的受益和风险曲线。如图所示。,对于问题（2），我们用在01内的200个等分点，求出最优投资中的受益和风险曲线。如图所示。,（3）投资者的性格特征,风险收益曲线是上升曲线，且是向上凸的,从上图可看出，风险收益曲线是离散的。问题（1）只有5个点，问题（2）只有13个点。,这说明，在变化时，共有数次跳变，使收益和风险发生突变。,以问题（1）为例说明,当在00.0375变化时，投资者比较侧重于不冒风险，故给出的解为全部存入银行。我们称这类投资者为极憎风险者。,当在0.03800.2330变化时，由于投资对风险的要求比较平衡，故尽管存在利润的微小差异，但总体上讲比较稳定，我们称这些投资者为风险折衷者。,当在0.23351变化时，由于投资十分注重利益的获得，所以在投资时，他将所有资金投入收益最高的项目，我们称这些投资者为风险投资者。,对于一般情况进行观察后，我们发现当0.2时，若发生微小变化，投资额将变得越来越分散，这一点可从数据表中看出。这与分散投资以避免风险的方法是相一致的。,对于一般情况，我们可以用同样的分析方法以求得突变点，从而将投资这分开。,模型的灵敏度分析,基于投资者的性格特征，本模型中存在几个跳变点，在跳变点附近，利润与风险会有大的变化，它主要的缘于投资者的心理变化，我们姑且给它们一个名字心理极限点，在两个极限点之间，当值作微小变化，我们可以发现投资额并无太大变化，所以我们认为本模型在投资额定下后比较稳定，即微小的心理变化不足以扰动投资额。,模型的改进与推广,（1）在不需要太高的精度而只需要大致比例时，在分配资金以前可根据qi与qj的比例，确定两者资金比例，即,这样可大大简化算法。这种方法称为等风险分散投资。,（2）公司往往会对收益和风险提出一定的要求，在本模型中只需增加约束即可。,