多目标规划培训教材-课件

MATLAB求解多目标规划江西师范大学数信学院吴根秀MATLAB求解多目标规划江西师范大学数信学院吴根秀一、0-1规划的MATLAB求解数学模型：MINf’xS.T.Ax<=bAeqx=beqx=0,1命令格式：x=bintprog(f)x=bintprog(f,A,b)x=bintprog(f,A,b,Aeq,beq)x=bintprog(f,A,b,Aeq,beq,x0)x=bintprog(f,A,b,Aeq,beq,x0,options)[x,fval]=bintprog(...)[x,fval,exitflag]=bintprog(...)[x,fval,exitflag,output]=bintprog(...)一、0-1规划的MATLAB求解数学模型：MINf’x数学模型：MINlambdaS.T.F(x)-weight*lambda<=goal(达到目标）Ax<=b（线性不等式约束）Aeqx=beq（线性等式约束）C(x)<=0（非线性不等式约束）Ceq(x)=0（非线性等式约束）lb<=x<=ubF=[f1(x),f2(x),…]为多目标的目标函数；F与[C(x),Ceq(x)]都是通过function来定义；命令格式：x=fgoalattain(fun,x0,goal,weight)x=fgoalattain(fun,x0,goal,weight,A,b)x=fgoalattain(fun,x0,goal,weight,A,b,Aeq,beq)x=fgoalattain(fun,x0,goal,weight,A,b,Aeq,beq,lb,ub)二、多目标规划的MATLAB求解数学模型：MINlambda二、多目标规划的MATLAB求命令格式：x=fgoalattain(fun,x0,goal,weight)x=fgoalattain(fun,x0,goal,weight,A,b)x=fgoalattain(fun,x0,goal,weight,A,b,Aeq,beq)x=fgoalattain(fun,x0,goal,weight,A,b,Aeq,beq,lb,ub)x=fgoalattain(fun,x0,goal,weight,A,b,Aeq,beq,lb,ub,nonlcon)x=fgoalattain(fun,x0,goal,weight,A,b,Aeq,beq,...lb,ub,nonlcon,options)[x,fval]=fgoalattain(...)[x,fval,attainfactor]=fgoalattain(...)[x,fval,attainfactor,exitflag]=fgoalattain(...)[x,fval,attainfactor,exitflag,output]=fgoalattain(...)[x,fval,attainfactor,exitflag,output,lambda]=fgoalattain(...)二、多目标规划的MATLAB求解命令格式：二、多目标规划的MATLAB求解x=fgoalattain(@myfun,x0,goal,weight,A,b,Aeq,beq,...lb,ub,@mycon)wheremyconisaMATLABfunctionsuchasfunction[c,ceq]=mycon(x)c=...%computenonlinearinequalitiesatx.ceq=...%computenonlinearequalitiesat二、多目标规划的MATLAB求解x=fgoalattain(@myfun,x0,goal,weight)wheremyfunisaMATLABfunctionsuchasfunctionF=myfun(x)F=...%Computefunctionvaluesatx.x=fgoalattain(@myfun,x0,goal有关优化参数设置：options=optimset(‘GradObj’,‘on’)目标函数的梯度方向参数设置为‘on’时，用下列函数定义：function[F,G]=myfun(x)F=...%Computethefunctionvaluesatxifnargout>1%TwooutputargumentsG=...%GradientsevaluatedatxEndThegradientconsistsofthepartialderivativedF/dxofeachFatthepointx.二、多目标规划的MATLAB求解有关优化参数设置：二、多目标规划的MATLAB求解二、多目标规划的MATLAB求解有关优化参数设置：options=optimset(‘GradConstr’,‘on’）约束条件的梯度方向参数设置为‘on’时，用下列函数定义：function[c,ceq,GC,GCeq]=mycon(x)c=...%Nonlinearinequalitiesatxceq=...%Nonlinearequalitiesatxifnargout>2%Nonlconcalledwith4outputsGC=...%GradientsoftheinequalitiesGCeq=...%GradientsoftheequalitiesEnd注意：一般weight=abs(goal)二、多目标规划的MATLAB求解有关优化参数设置：模型：x’=(A+BKC)x+Bu，设计K满足目标：Y=Cx1）循环系统的特征值（由命令eig(A+B*K*C)确定）的目标为goal=[-5,-3,-1]2）K中元素均在[-4,4]中;设特征值的weight=abs(goal),定义目标函数F如下：functionF=eigfun(K,A,B,C)F=sort(eig(A+B*K*C));%Evaluateobjectives，由小到大排列优化程序为：A=[-0.500;0-210;01-2];B=[10;-22;01];C=[100;001];K0=[-1-1;-1-1];%Initializecontrollermatrixgoal=[-5-3-1];%Setgoalvaluesfortheeigenvaluesweight=abs(goal)%Setweightforsamepercentagelb=-4*ones(size(K0));%Setlowerboundsonthecontrollerub=4*ones(size(K0));%Setupperboundsonthecontrolleroptions=optimset('Display','iter');%Setdisplayparameter[K,fval,attainfactor]=fgoalattain(@(K)eigfun(K,A,B,C)...goal,weight,[],[],[],[],lb,ub,[],options)二、举例---有关循环控制系统优化问题模型：x’=(A+BKC)x+Bu，设计K满足目标：二、举例运行结果如下Activeconstraints:124910K=-4.0000-0.2564-4.0000-4.0000fval=-6.9313-4.1588-1.4099attainfactor=-0.3863二、举例---有关循环控制系统优化问题运行结果如下二、举例---有关循环控制系统优化问题如果至少保证38.63%的目标精确匹配，设置‘GoalsExactAchieve’参数值为3options=optimset('GoalsExactAchieve',3);[K,fval,attainfactor]=fgoalattain(...@(K)eigfun(K,A,B,C),K0,goal,weight,[],[],[],[],lb,ub,[],...options)Afteraboutseveniterations,asolutionisK=-1.59541.2040-0.4201-2.9046fval=-5.0000-3.0000-1.0000attainfactor=1.0859e-20表明目标已完全匹配二、举例---有关循环控制系统优化问题如果至少保证38.63%的目标精确匹配，设置‘GoalsEx谢谢！谢谢！

初等模型举例多目标规划培训教材-课件

常见类型定性模型经验公式（拟合、插值）量纲分析比例模型常见类型定性模型§2.1

崖高的估算假如你站在崖顶且身上带着一只具有跑表功能的计算器，你也许会出于好奇心想用扔下一块石头听回声的方法来估计山崖的高度，假定你能准确地测定时间，你又怎样来推算山崖的高度呢，请你分析一下这一问题。我有一只具有跑表功能的计算器。§2.1崖高的估算假如你站在崖顶且身上方法一假定空气阻力不计，可以直接利用自由落体运动的公式来计算。例如，设t=4秒，g=9.81米/秒2，则可求得h≈78.5米。

我学过微积分，我可以做得更好，呵呵。

方法一假定空气阻力不计，可以直接利用自由落体运动的公式我学除去地球吸引力外，对石块下落影响最大的当属空气阻力。根据流体力学知识，此时可设空气阻力正比于石块下落的速度，阻力系数K为常数，因而，由牛顿第二定律可得：

令k=K/m，解得

代入初始条件v(0)=0，得c=－g/k，故有

再积分一次，得：

除去地球吸引力外，对石块下落影响最大的当属空气阻力。根若设k=0.05并仍设t=4秒，则可求得h≈73.6米。

听到回声再按跑表，计算得到的时间中包含了反应时间

进一步深入考虑不妨设平均反应时间为0.1秒，假如仍设t=4秒，扣除反应时间后应为3.9秒，代入式①，求得h≈69.9米。

①多测几次，取平均值再一步深入考虑代入初始条件h(0)=0，得到计算山崖高度的公式：

将e-kt用泰勒公式展开并令k→0+

，即可得出前面不考虑空气阻力时的结果。若设k=0.05并仍设t=4秒，则可求得h≈73.6米。还应考虑回声传回来所需要的时间。为此，令石块下落的真正时间为t1，声音传回来的时间记为t2，还得解一个方程组：这一方程组是非线性的，求解不太容易，为了估算崖高竟要去解一个非线性主程组似乎不合情理

相对于石块速度，声音速度要快得多，我们可用方法二先求一次

h，令t2=h/340，校正t，求石块下落时间t1≈t-t2将t1代入式①再算一次，得出崖高的近似值。例如，若h=69.9米，则t2≈0.21秒，故t1≈3.69秒，求得h≈62.3米。还应考虑回声传回来所需要的时间。为此，令石块下落的真§2.2

录像带还能录多长时间录像机上有一个四位计数器，一盘180分钟的录像带在开始计数时为0000，到结束时计数为1849，实际走时为185分20秒。我们从0084观察到0147共用时间3分21秒。若录像机目前的计数为1428，问是否还能录下一个60分钟的节目？§2.2录像带还能录多长时间录像机上有一个四位计数器，一rθRl由得到又因和得

积分得到即从而有我们希望建立一个录像带已录像时间t与计数器计数n之间的函数关系。为建立一个正确的模型，首先必须搞清哪些量是常量，哪些量是变量。首先，录像带的厚度W是常量，它被绕在一个半径为r的园盘上，见图。磁带转动中线速度v显然也是常数，否则图象声音必然会失真。此外，计数器的读数n与转过的圈数有关，从而与转过的角度θ成正比。rθRl由得到又rθRl

此式中的三个参数W、v和r均不易精确测得，虽然我们可以从上式解出t与n的函数关系，但效果不佳，故令则可将上式简化为：故令上式又可化简记成t=an2+bn

rθRl此式中的三个参数W、v和r均不易精确测得，虽然我们t=an2+bn

rθRl上式以a、b为参数显然是一个十分明智的做法，它为公式的最终确立即参数求解提供了方便。将已知条件代入，得方程组：从后两式中消去t1，解得a=0.0000291,b=0.04646,故t=0.0000291n2+0.04646n，令n=1428，得到t=125.69（分）由于一盒录像带实际可录像时间为185.33分，故尚可录像时间为59.64分，已不能再录下一个60分钟的节目了。

t=an2+bnrθRl上式以a、b为参数显然是一个十分在解决实际问题时，注意观察和善于想象是十分重要的，观察与想象不仅能发现问题隐含的某些属性，有时还能顺理成章地找到解决实际问题的钥匙。本节的几个例子说明，猜测也是一种想象力。没有合理而又大胆的猜测，很难做出具有创新性的结果。开普勒的三大定律（尤其是后两条）并非一眼就能看出的，它们隐含在行星运动的轨迹之中，隐含在第谷记录下来的一大堆数据之中。历史上这样的例子实在太多了。在获得了一定数量的资料数据后，人们常常会先去猜测某些结果，然后试图去证明它。猜测一经证明就成了定理，而定理一旦插上想象的翅膀，又常常会被推广出许多更为广泛的结果。即使猜测被证明是错误的，结果也决不是一无所获的失败而常常是对问题的更为深入的了解。§2.3最短路径与最速方案问题在解决实际问题时，注意观察和善于想象是十分重要的，观察与想象例5（最短路径问题）设有一个半径为r的圆形湖，圆心为O。A、B

位于湖的两侧，AB连线过O，见图。现拟从A点步行到B点，在不得进入湖中的限制下，问怎样的路径最近。

ABOr将湖想象成凸出地面的木桩，在AB间拉一根软线，当线被拉紧时将得到最短路径。根据这样的想象，猜测可以如下得到最短路径：过A作圆的切线切圆于E，过B作圆的切线切圆于F。最短路径为由线段AE、弧EF和线段FB连接而成的连续曲线（根据对称性，AE′，弧E′F′，F′B连接而成的连续曲线也是）。EFE′F′例5（最短路径问题）设有一个半径为r的圆形湖，圆心为以上只是一种猜测，现在来证明这一猜测是正确的。为此，先介绍一下凸集与凸集的性质。定义2.1（凸集）称集合R为凸集，若x1、x2∈R及λ∈[0，1]，总有λx1+（1+λ）x2∈R。即若x1、x2∈R，则x1、x2的连线必整个地落在R中。定理2.2（分离定理）对平面中的凸集R与R外的一点K，存在直线l,l

分离R与K，即R与K分别位于l的两侧（注：对一般的凸集R与R外的一点K，则存在超平面分离R与K），见图。klR下面证明猜想以上只是一种猜测，现在来证明这一猜测是正确的。为此，先介绍一猜测证明如下：（方法一）显然，由AE、EF、FB及AE′，E′F′，F′B围成的区域R是一凸集。利用分离定理易证最短径不可能经过R外的点，若不然，设Γ为最短路径，Γ过R外的一点M，则必存在直线l分离M与R，由于路径Γ是连续曲线，由A沿Γ到M，必交l于M1，由M沿Γ到B又必交l于M2。这样，直线段M1M2的长度必小于路径M1MM2的长度，与Γ是A到B的最短路径矛盾，至此，我们已证明最短路径必在凸集R内。不妨设路径经湖的上方到达B点，则弧EF必在路径F上，又直线段AE是由A至E的最短路径，直线FB是由F到B的最短路径，猜测得证。ABOrEFE′F′M1M2MΓl猜测证明如下：（方法一）显然，由AE、EF、FB及AE′，还可用微积分方法求弧长，根据计算证明满足限止条件的其他连续曲线必具有更大的长度；此外，本猜测也可用平面几何知识加以证明等。根据猜测不难看出，例5中的条件可以大大放松，可以不必设AB过圆心，甚至可不必设湖是圆形的。例如对下图，我们可断定由A至B的最短路径必为l1与l2之一，其证明也不难类似给出。ABl1l2D到此为止，我们的研讨还只局限于平面之中，其实上述猜测可十分自然地推广到一般空间中去。1973年，J.W.Craggs证明了以上结果：若可行区域的边界是光滑曲面。则最短路径必由下列弧组成，它们或者是空间中的自然最短曲线，或者是可行区域的边界弧。而且，组成最短路径的各段弧在连接点处必定相切。还可用微积分方法求弧长，根据计算证明满足限止条件的其他连续曲例6

一辆汽车停于A处并垂直于AB方向，此汽车可转的最小圆半径为R，求不倒车而由A到B的最短路径。解（情况1）若|AB|>2R，最短路径由弧AC与切线BC组成（见图①

）。（情况2）若|AB|<2R，则最短路径必居于图②（a）、（b）两曲线之中。可以证明，（b）中的曲线ABC更短。AR2RBRC①②ABoC（a）CABo1o2（b）例6一辆汽车停于A处并垂直于AB方向，此解（情况1）若例7

驾驶一辆停于A处与AB成θ1角度的汽车到B处去，已知B处要求的停车方向必须与AB成θ2角，试找出最短路径（除可转的最小圆半径为R外，不受其他限止）。解根据Craggs定理并稍加分析可知，最短路径应在l1与l2中，见图，比较l1与l2的长度，即可得到最短路径。Al1l2Bθ2θ1例7驾驶一辆停于A处与AB成θ1角度的汽解根据Cra最速方案问题例8

将一辆急待修理的汽车由静止开始沿一直线方向推至相隔S米的修车处，设阻力不计，推车人能使车得到的推力f满足:-B≤f≤A,f>0为推力，f<0为拉力。问怎样推车可使车最快停于修车处。

设该车的运动速度为υ=υ(t)，根据题意，υ(0)=υ(T)=0，其中T为推车所花的全部时间。由于-B≤f≤A，且f=mυ′，可知-b≤υ′≤a（其中m为汽车质量，a=A/m,b=B/m）。据此不难将本例归纳为如下的数学模型：

minT

υ(0)=υ(T)=0最速方案问题例8将一辆急待修理的汽车由静止开始沿一此问题为一泛函极值问题，求解十分困难，为得出一个最速方案。我们作如下猜测：猜测最速方案为以最大推力将车推到某处，然后以最大拉力拉之，使之恰好停于修车处，其中转换点应计算求出证明设υ=υ(t)为在最速推车方案下汽车的速度，则有。设此方案不同于我们的猜测。现从O点出发，作射线y=at；从（t,0）出发，作直线y=-b(t-T)交y=at于A，由于，曲线υ=υ(t)必位于三角形区域DAT的内部，从而有ΔOAT的面积大于S。在O到T之间任取一点T′

，过T′作AT的平行线交OA于A′

。显然ΔOA′T′的面积S（T′）是T′的连续函数，当T′=0时S（0）=0，当T′=T时，S（T）>S，故由连续函数的性质存在某T′<T，S（T′）=S但这一结果与υ=υ(t)是最优方案下的车速的假设矛盾，因为用我们猜测的推车方法推车，只需T′时间即可将车推到修车处，而T′<T。oυtAT′TA′Sy=aty=-b(t-T)此问题为一泛函极值问题，求解十分困难，为得出一个最速方案。我逻辑模型多目标规划培训教材-课件例1

拟将一批尺寸为1×2×4的的商品装入尺寸为6×6×6的正方体包装箱中，问是否存在一种装法，使装入的该商品正好充满包装箱。解

将正方体剖分成27个2×2×2的小正方体，并按下图所示黑白相间地染色。再将每一2×2×2的小正方体剖分成1×1×1的小正方体。易见，27个2×2×2的正方体中，有14个是黑的，13个是白的（或13黑14白），故经两次剖分，共计有112个1×1×1的黑色小正方体和104个1×1×1的白色小正方体。虽然包装箱的体积恰好是商品体积的27倍，但容易看到，不论将商品放置在何处，它都将占据4个黑色和4个白色的1×1×1小正方体的位置，故商品不可能充满包装箱。例1拟将一批尺寸为1×2×4的的商品装入尺寸为6×6

德国著名的艺术家AlbrechtDürer(1471-1521)于1514年曾铸造了一枚名为“MelencotiaI”的铜币。令人奇怪的是在这枚铜币的画面上充满了数学符号、数字及几何图形。这里，我们仅研究铜币右上角的数字问题

例2.Dürer魔方（或幻方）问题德国著名的艺术家AlbrechtDürer(147所谓的魔方是指由1~n2这n2个正整数按一定规则排列成的一个n行n列的正方形。n称为此魔方的阶。Dürer魔方：4阶，每一行之和为34，每一列之和为34，对角线（或反对角线）之和是34，每个小方块中的数字之和是34，四个角上的数字加起来也是34什么是Dürer魔方多么奇妙的魔方！铜币铸造时间：1514年

所谓的魔方是指由1~n2这n2个正整数按一定规则排列成的一个

构造魔方是一个古老的数学游戏，起初它还和神灵联系在一起，带有深厚的迷信色彩。传说三千二百多年前（公元前2200年），因治水出名皇帝大禹就构造了三阶魔方（被人们称“洛书”），至今还有人把它当作符咒用于某些迷信活动，大约在十五世纪时，魔方传到了西方，著名的科尼利厄斯·阿格里帕（1486-1535）先后构造出了3~9阶的魔方。构造魔方是一个古老的数学游戏，起初它还和神灵联系在如何构造魔方奇数（不妨n=5）阶的情况Step1:在第一行中间写1Step2:每次向右上方移一格依次填按由小到大排列的下一个数，向上移出界时填下一列最后一行的小方格；向右移出界时填第一列上一行的小方格。若下面想填的格已填过数或已达到魔方的右上角时，改填刚才填的格子正下方的小方格，继续Step2直到填完12345678910111213141516171819202122232425偶数阶的情况

偶数阶的魔方可以利用奇数阶魔方拼接而成，拉尔夫·斯特雷奇给出了一种拼接的方法，这里不作详细介绍如何构造魔方奇数（不妨n=5）阶的情况Step1:在第五阶没人知道有多少个！！！三阶1个反射和中心旋转生成8个四阶880个反射和中心旋转生成7040个魔方数量随阶数n增长的速度实在是太惊人了！同阶魔方的个数五阶没人知道有多少个！！！三阶允许构成魔方的数取任意实数允许取实数

n阶魔方A、B，任意实数α、βαA+βB是n阶魔方具有指定性质的魔方全体构成一个线性空间问题已发生了实质性变化注：刻画一个线性空间只需指出它的维数并求出此线性空间的一组基底松驰问题的讨论允许构成魔方的数取任意实数允许取实数n阶魔方A、

1在第一行中共有4种取法，为保持上述性质的成立，第二行中的1还有两种取法。当第二行的1也取定后，第三行与第四行的1就完全定位了，故一共可作出8个不同的最简方阵，称之为基本魔方并记之为Q1，…，Q8

仍以4阶方阵为例。令R为行和，C为列和，D为对角线和，S为小方块和定义0-方：R=C=D=S=0定义1-方：R=C=D=S=4R=C=D=S=1的方阵构成的线性空间具有什么样的性质？类似于构造n维欧氏空间的标准基，利用0和1我们来构造一些R=C=D=S=1的最简单的方阵。1在第一行中共有4种取法，为保持上述性质的

显然，Dürer空间（简称D空间）中任何一个元素都可以用Q1，Q2，…，Q8来线性表示，但它们能否构成D空间的一组基呢？

显然，Dürer空间（简称D空间）中任何一容易看出：Q1，…，Q8这8个基本方是线性相关的，即至少存在一个Qj，可以通过其它7个基本方的线性组合得到。这8个基本方的地位是等同的，故可不妨设j=8。下面验证Q1，Q2，…，Q7是否线性相关。

令：，即=容易看出：Q1，…，Q8这8个基本方是线性相关的，即至少存在等号两边对应元素相比较，得r1=r2=…=r7=0，所以是线性无关是D空间的最小生成集。

令D

即解方程组：

=

解得

D=研究AlbrechtDürer铸造的铜币等号两边对应元素相比较，得r1=r2=…=r7=0，是D空间2004年浙江大学数学建模竞赛（B题）通讯卫星上的开关设置

地面上存在着n个接收站与n个发送站，而在通讯卫星上则设置了若干种开关模式。开关模式可用矩阵P=(pij)来表示，若卫星可接收发送站i发射的信息并将信息传送回地面的接收站j时，矩阵中的元素pij=1，否则pij=0。通讯卫星上的接收发送任务也可以用一个矩阵T=（tij）来表示，其元素tij为需经通讯卫星传递的由i发点发送到j接收点的信息量的传送时间长度。由于技术上的原因，当发送站i在发送给接收站j信息时，它不能同时发送给别的接收站信息；同样，当接收站j在接收发送站i的信息时，也不能同时接收其他发送站发送的信息。你的任务是：2004年浙江大学数学建模竞赛

设计一组开关模式，k=1,…,r（注：r应当尽可能小），使得对任意给定的任务矩阵T，卫星开关设置均能完成要求的发接收任务。设计一个算法，在发接收任务T给出后，可根据你设计的开关模式（k=1,…,r）求出各开关的使用时间λk，使得在完成预定传送任务的前提下使用各开关模式的总时间最短。同样由于技术上的原因，开关模式的总数r有一个上限。当需要传送的任务数数量较大时，仍无法分派任务。请你想一些办法来解决这一困难，（当然，这时你可能要作出一些牺牲，即传送时间可能会增加一些）。

设计一组开关模式，k=1,…,r（注：r应当尽可能小）问题及模型问题的标准形式为：在地面上存在着n个收站与n个发战，而在通讯卫星上则设置了若干种开关模式。开关模式可用矩阵P=(pij)来表示，若卫星可接收发射站i发射的信息并将信息传送回地面的接收站j时，矩阵元素pij=1，否则pij=0。通讯卫星的接发任务也可用一矩阵T=（tij）来表示，其元素tij为需经通讯卫星传递的由i发点发送到j接受点的信息量的传送时间长度。问题要求求r并设计一组开关模式Pk，k=1,…,r及模式Pk的使用时间λk，使得在完成预定传送任务的前提下各开关模式使用的总时间最短，即要求求解下面的问题：问题及模型问题的标准形式为：在地面上存在着n个收站与n个发战例1

设

这是一个有3个发送站与3个接收站的实例，tij在矩阵中已给出，例如由发站1传送到收站1的通讯量为3单位时间等。分析

容易看出，三个发站需传送的时间分别为6、5、5；而三个收站需接收的时间分别为6、3、7。为完成全部传送任务，通讯卫星总传送时间至少应为7单位时间，即的下界为7。由于技术上的原因，当发站i在发送给收站j信息时，它不能同时发送给别的收站信息；同样，当收站j在接收发站i的信息时，也不能同时接收其他发站发送的信息。这一要求说明，任一开关模式Pk应具有以下性质：（1）Pk的每一行中有且只有一个1，每一列中也有且只有一个1；（2）所有的1均位于不同的行列中。满足（1）、（2）的矩阵被称为置换矩阵，n阶置换矩阵Pk共有n!个，当n较大时，我们不可能在通讯卫星上设置这么多种不同的开关模式。因而，为了设计出切实可行的开关模式，我们还得另想办法。例1设这是一个有3个发送站与3个接收站的实例，tij在（设计方法1）注意到Pk每行（或列）元素之和均为1，故不管如何指派开关的使用时间（即不论如何取λk），矩阵均具有某些特殊的性质，例如其行和（及列和）均为同一常数。这样的矩阵构成一个线性空间（参见逻辑模型第一节Dürer魔方），为减少开关模式的种类，可取此空间的一组基底作为开关模式。在使用这种开关模式时，无论T的元素tij怎么取，通讯卫星对每一发（收）点的开通时间总和是恒定的。在这种开关模式下，可按如下方式指派各开关模式的使用时间：步1先将T改变为，满足：（1）≥T（2）记，步2用Pk表示，即将分解为（r为空间的维数）（设计方法1）注意到Pk每行（或列）元素之和均为1，故不管如将T化为的方法一般有无穷多种，如可如下化法：令事实上，，（即通讯卫星传送总时间的下界）。令其中用这种方法化例中的T，得到的任一行（或列）中元素之和均为7。将T化为的方法一般有无穷多种，如可如下化法：令定义1称行和、列%%和均相等的矩阵为双随机矩阵（Doublystochasticmatrix）定理1

（Birkhoff定理，1944）任一n阶双随机矩阵均可写成至多（n－1）2+1个置换矩阵的非线性组合。的分解可如下进行：步1选取由Pij>0可推出>0的置换矩阵P步2确定

步3取，用－代替步4若=0，停；否则，返回步1。例2.

为方便起见，我们来分解一个元素均为非负整数的3阶双随机矩阵，（由Birkhoff定理，r≤5）定义1定理1（Birkhoff定理，1944）任一n阶双随解：取，λ=min{1,3,3}=1分解成，再取

因min{5,5,3}=3，又有，取

解：取于是又有易得分解结果为：于是又有易得分解结果为：尚需解决的问题是如何求P，使得Pij>0必有。读者不难发现，此问题可以通过求解一个两分图上的最大流（或最大匹配）来实现，计算量为O(n4)，是多项式时间可解的。具体方法为：作一两分图，若，则作边（i,j），令边容量为1，这样，可作出P的充要条件是该最大流问题的最大流量为n。对例9.33，n=3。由于所有，先取

，－P1为

相应的两分图为：于是又可求得尚需解决的问题是如何求P，使得Pij>0必有，相应的两分图为：又可得，…，如此下去，直到作不出P为至，由于的特殊性质及Birkhoff定理，上述分解必能在不超过r=(n－1)2+1步内终止。上述开关设计方法要求在通讯卫星上设置(n－1)2+1种不同的开关模式（即Pk），当n稍大时，(n－1)2+1仍显得太大而使得使用时不便。例如，当n=41时，(n－1)2+1=|60|。为实用方便，人们研究了限止开关模式个数的相应问题。，相应的两分图为：又可得若要求r≤n，即要求通讯卫星上至多设置n种开关模式，则问题化为令r≤n，求不超过n个置换矩阵Pk及λk，使之满足：minS.t

为了使任意一对发射法与接收站之间的传送均为可能实现的，自然应要求Pk满足（5.1）

（5.2）

（右面的矩阵有n2个值为1的分量，每一Pk

恰有n个1分量）故r=n。容易看出，（5.1）隐含着T的每一元素只能被唯一的P复盖，即T的元素在分解中是不可分割的，这当然是一个好性质，使实际操作时较为方便，但可惜的是对一般的双随机矩阵，分解很可能无解。若要求r≤n，即要求通讯卫星上至多设置n种开关模式，则问题化例3

若取，

，

（注意：T已是双随机矩阵，行和列和均为10）则min

S.t

的解为λ1=3，λ2=4，λ3=5。例3若取，，（注意：T已是双随机矩阵，行和列和均为1（大于10）而但等号经常并不成立。1985年，F·Rendel证明，在给定满足（5.2）的置换矩阵P1,…,Pn后，求解问题（5.1）是NP难的，从而不可能存在多项式时间算法，除非P=NP。现要求r≤2n一种自然而方便的开关设置为引入两组各有n个开关模式的置换矩阵P1,…,Pn，Q1,…,Qn，满足下面的（5.3)式：例如，当n=3时，可令：（大于10）而但等号经常并不成立。1985年，F·Rend（注：这种设置方法保持了其内在的对称性，不失为一种明智的做法。）现在，我们来分解例9.33中的双随机矩阵，令=，得方程组（注：这种设置方法保持了其内在的对称性，不失为一种明智的做法求出各对角线与反对角线上的三个元素之和，并作一些简单的消去运算；将矩阵的所有元素相加，可得下面的方程组：注意到（5.3），易证空间的维数为5，故之一可任取，（稍加注意即可保持非负性），例如，令μ3=0，求得，故有求出各对角线与反对角线上的三个元素之和，并作一些简单的消去运读者不难验证，上述方法可推广到n是奇数的一般情况。事实，由各对角线元素之和可导出n－1个方程，由各反对角线元素之和又可导出n－1个方程，加上矩阵所有元素之和导出的等式，共计可导出2n－1个方程，并易知它们是独立的。另一方面空间的维数恰为2n－1，故之一可任取，而通过方程组解得所有的，（只须注意保持其非负性即可）但当n为偶数时，情况就不大相同了。让我们先来观察一下n=4的情况。当n=4时，读者不难验证，上述方法可推广到n是奇数的一般情况。事实，由各易见，具有非常特殊的结构，一般的偶数阶双随机矩阵，即使其元素是非负整数，也无法用Pk、Qk来分解。当具有上述结构时，能否用Pk和Qk来分解呢？易见，由各对角线元素之和可导出：易见，另外，由反对角线元素之和又可导出上述方程中只有6个是独立的，且已不可能再得出新的独立方程，（读者可自行分析之）故可选取其中2个的值，进而可解出其余。例如，若令

4=3=0，可得2=1，1=0，进而可求得1=2，4=3，3=3及2=4，

已达到下界。易见，P1+P3=Q2+Q4，P2+P4=Q1+Q3，故空间的维数为6，与上面的分析是一致的。读者可将上述讨论推广到n为一般偶数的情况，分析方法是完全类似的。另外，由反对角线元素之和又可导出上述方程中只有6个是独立的，当n是偶数时，我们虽无法将一般的双随机矩阵分解为Pk

、Qk的非负组合，但上述讨论仍然是十分有意义的。首先，要求完成的任务矩阵是T，在将T转换成时我们可尽量使具有上述的特殊结构（有兴趣的读者可自行研究这一问题），只要能做到这一点，即可给出一个达到下界的开关模式的指派方式。其次，即使这样的努力没有成功，也容易给出一个具有上述特殊结构矩阵，并使尽可能地小，即给出一种开关指派的近似最佳方法，由此可设计出效果较好的近似算法。由于技术水平的提高，目前通讯卫星传送信息已允许一个发射站同时向多个接收站发送信息，当然，同时发送的信息条数具有某一上限，例如上限为v。1987年，J.L.Lewandowski和C.L.Liu研究了如下更一般的问题：当n是偶数时，我们虽无法将一般的双随机矩阵分解为Pk、Qk给定一正整数v，（v为通讯卫星传送容量的总限止），求开关模式M：=：={;(0,1)|,i=1,…,m;,i=1,…,n,}的设计，要求所用的开关模式总数量r尽可能小，并使有解，其中T为信息传送量矩阵（需满足一定要求），ak为开关模式Mk的使用时间。他们设计了一个求解此问题的O(n5)算法，有兴趣的读者可直接阅读他们的论文。

TheEnd给定一正整数v，（v为通讯卫星传送容量的总限止），求开关模式MATLAB求解多目标规划江西师范大学数信学院吴根秀MATLAB求解多目标规划江西师范大学数信学院吴根秀一、0-1规划的MATLAB求解数学模型：MINf’xS.T.Ax<=bAeqx=beqx=0,1命令格式：x=bintprog(f)x=bintprog(f,A,b)x=bintprog(f,A,b,Aeq,beq)x=bintprog(f,A,b,Aeq,beq,x0)x=bintprog(f,A,b,Aeq,beq,x0,options)[x,fval]=bintprog(...)[x,fval,exitflag]=bintprog(...)[x,fval,exitflag,output]=bintprog(...)一、0-1规划的MATLAB求解数学模型：MINf’x数学模型：MINlambdaS.T.F(x)-weight*lambda<=goal(达到目标）Ax<=b（线性不等式约束）Aeqx=beq（线性等式约束）C(x)<=0（非线性不等式约束）Ceq(x)=0（非线性等式约束）lb<=x<=ubF=[f1(x),f2(x),…]为多目标的目标函数；F与[C(x),Ceq(x)]都是通过function来定义；命令格式：x=fgoalattain(fun,x0,goal,weight)x=fgoalattain(fun,x0,goal,weight,A,b)x=fgoalattain(fun,x0,goal,weight,A,b,Aeq,beq)x=fgoalattain(fun,x0,goal,weight,A,b,Aeq,beq,lb,ub)二、多目标规划的MATLAB求解数学模型：MINlambda二、多目标规划的MATLAB求命令格式：x=fgoalattain(fun,x0,goal,weight)x=fgoalattain(fun,x0,goal,weight,A,b)x=fgoalattain(fun,x0,goal,weight,A,b,Aeq,beq)x=fgoalattain(fun,x0,goal,weight,A,b,Aeq,beq,lb,ub)x=fgoalattain(fun,x0,goal,weight,A,b,Aeq,beq,lb,ub,nonlcon)x=fgoalattain(fun,x0,goal,weight,A,b,Aeq,beq,...lb,ub,nonlcon,options)[x,fval]=fgoalattain(...)[x,fval,attainfactor]=fgoalattain(...)[x,fval,attainfactor,exitflag]=fgoalattain(...)[x,fval,attainfactor,exitflag,output]=fgoalattain(...)[x,fval,attainfactor,exitflag,output,lambda]=fgoalattain(...)二、多目标规划的MATLAB求解命令格式：二、多目标规划的MATLAB求解x=fgoalattain(@myfun,x0,goal,weight,A,b,Aeq,beq,...lb,ub,@mycon)wheremyconisaMATLABfunctionsuchasfunction[c,ceq]=mycon(x)c=...%computenonlinearinequalitiesatx.ceq=...%computenonlinearequalitiesat二、多目标规划的MATLAB求解x=fgoalattain(@myfun,x0,goal,weight)wheremyfunisaMATLABfunctionsuchasfunctionF=myfun(x)F=...%Computefunctionvaluesatx.x=fgoalattain(@myfun,x0,goal有关优化参数设置：options=optimset(‘GradObj’,‘on’)目标函数的梯度方向参数设置为‘on’时，用下列函数定义：function[F,G]=myfun(x)F=...%Computethefunctionvaluesatxifnargout>1%TwooutputargumentsG=...%GradientsevaluatedatxEndThegradientconsistsofthepartialderivativedF/dxofeachFatthepointx.二、多目标规划的MATLAB求解有关优化参数设置：二、多目标规划的MATLAB求解二、多目标规划的MATLAB求解有关优化参数设置：options=optimset(‘GradConstr’,‘on’）约束条件的梯度方向参数设置为‘on’时，用下列函数定义：function[c,ceq,GC,GCeq]=mycon(x)c=...%Nonlinearinequalitiesatxceq=...%Nonlinearequalitiesatxifnargout>2%Nonlconcalledwith4outputsGC=...%GradientsoftheinequalitiesGCeq=...%GradientsoftheequalitiesEnd注意：一般weight=abs(goal)二、多目标规划的MATLAB求解有关优化参数设置：模型：x’=(A+BKC)x+Bu，设计K满足目标：Y=Cx1）循环系统的特征值（由命令eig(A+B*K*C)确定）的目标为goal=[-5,-3,-1]2）K中元素均在[-4,4]中;设特征值的weight=abs(goal),定义目标函数F如下：functionF=eigfun(K,A,B,C)F=sort(eig(A+B*K*C));%Evaluateobjectives，由小到大排列优化程序为：A=[-0.500;0-210;01-2];B=[10;-22;01];C=[100;001];K0=[-1-1;-1-1];%Initializecontrollermatrixgoal=[-5-3-1];%Setgoalvaluesfortheeigenvaluesweight=abs(goal)%Setweightforsamepercentagelb=-4*ones(size(K0));%Setlowerboundsonthecontrollerub=4*ones(size(K0));%Setupperboundsonthecontrolleroptions=optimset('Display','iter');%Setdisplayparameter[K,fval,attainfactor]=fgoalattain(@(K)eigfun(K,A,B,C)...goal,weight,[],[],[],[],lb,ub,[],options)二、举例---有关循环控制系统优化问题模型：x’=(A+BKC)x+Bu，设计K满足目标：二、举例运行结果如下Activeconstraints:124910K=-4.0000-0.2564-4.0000-4.0000fval=-6.9313-4.1588-1.4099attainfactor=-0.3863二、举例---有关循环控制系统优化问题运行结果如下二、举例---有关循环控制系统优化问题如果至少保证38.63%的目标精确匹配，设置‘GoalsExactAchieve’参数值为3options=optimset('GoalsExactAchieve',3);[K,fval,attainfactor]=fgoalattain(...@(K)eigfun(K,A,B,C),K0,goal,weight,[],[],[],[],lb,ub,[],...options)Afteraboutseveniterations,asolutionisK=-1.59541.2040-0.4201-2.9046fval=-5.0000-3.0000-1.0000attainfactor=1.0859e-20表明目标已完全匹配二、举例---有关循环控制系统优化问题如果至少保证38.63%的目标精确匹配，设置‘GoalsEx谢谢！谢谢！

初等模型举例多目标规划培训教材-课件

常见类型定性模型经验公式（拟合、插值）量纲分析比例模型常见类型定性模型§2.1

崖高的估算假如你站在崖顶且身上带着一只具有跑表功能的计算器，你也许会出于好奇心想用扔下一块石头听回声的方法来估计山崖的高度，假定你能准确地测定时间，你又怎样来推算山崖的高度呢，请你分析一下这一问题。我有一只具有跑表功能的计算器。§2.1崖高的估算假如你站在崖顶且身上方法一假定空气阻力不计，可以直接利用自由落体运动的公式来计算。例如，设t=4秒，g=9.81米/秒2，则可求得h≈78.5米。

我学过微积分，我可以做得更好，呵呵。

方法一假定空气阻力不计，可以直接利用自由落体运动的公式我学除去地球吸引力外，对石块下落影响最大的当属空气阻力。根据流体力学知识，此时可设空气阻力正比于石块下落的速度，阻力系数K为常数，因而，由牛顿第二定律可得：

令k=K/m，解得

代入初始条件v(0)=0，得c=－g/k，故有

再积分一次，得：

除去地球吸引力外，对石块下落影响最大的当属空气阻力。根若设k=0.05并仍设t=4秒，则可求得h≈73.6米。

听到回声再按跑表，计算得到的时间中包含了反应时间

进一步深入考虑不妨设平均反应时间为0.1秒，假如仍设t=4秒，扣除反应时间后应为3.9秒，代入式①，求得h≈69.9米。

①多测几次，取平均值再一步深入考虑代入初始条件h(0)=0，得到计算山崖高度的公式：

将e-kt用泰勒公式展开并令k→0+

，即可得出前面不考虑空气阻力时的结果。若设k=0.05并仍设t=4秒，则可求得h≈73.6米。还应考虑回声传回来所需要的时间。为此，令石块下落的真正时间为t1，声音传回来的时间记为t2，还得解一个方程组：这一方程组是非线性的，求解不太容易，为了估算崖高竟要去解一个非线性主程组似乎不合情理

相对于石块速度，声音速度要快得多，我们可用方法二先求一次

h，令t2=h/340，校正t，求石块下落时间t1≈t-t2将t1代入式①再算一次，得出崖高的近似值。例如，若h=69.9米，则t2≈0.21秒，故t1≈3.69秒，求得h≈62.3米。还应考虑回声传回来所需要的时间。为此，令石块下落的真§2.2

录像带还能录多长时间录像机上有一个四位计数器，一盘180分钟的录像带在开始计数时为0000，到结束时计数为1849，实际走时为185分20秒。我们从0084观察到0147共用时间3分21秒。若录像机目前的计数为1428，问是否还能录下一个60分钟的节目？§2.2录像带还能录多长时间录像机上有一个四位计数器，一rθRl由得到又因和得

积分得到即从而有我们希望建立一个录像带已录像时间t与计数器计数n之间的函数关系。为建立一个正确的模型，首先必须搞清哪些量是常量，哪些量是变量。首先，录像带的厚度W是常量，它被绕在一个半径为r的园盘上，见图。磁带转动中线速度v显然也是常数，否则图象声音必然会失真。此外，计数器的读数n与转过的圈数有关，从而与转过的角度θ成正比。rθRl由得到又rθRl

此式中的三个参数W、v和r均不易精确测得，虽然我们可以从上式解出t与n的函数关系，但效果不佳，故令则可将上式简化为：故令上式又可化简记成t=an2+bn

rθRl此式中的三个参数W、v和r均不易精确测得，虽然我们t=an2+bn

rθRl上式以a、b为参数显然是一个十分明智的做法，它为公式的最终确立即参数求解提供了方便。将已知条件代入，得方程组：从后两式中消去t1，解得a=0.0000291,b=0.04646,故t=0.0000291n2+0.04646n，令n=1428，得到t=125.69（分）由于一盒录像带实际可录像时间为185.33分，故尚可录像时间为59.64分，已不能再录下一个60分钟的节目了。

t=an2+bnrθRl上式以a、b为参数显然是一个十分在解决实际问题时，注意观察和善于想象是十分重要的，观察与想象不仅能发现问题隐含的某些属性，有时还能顺理成章地找到解决实际问题的钥匙。本节的几个例子说明，猜测也是一种想象力。没有合理而又大胆的猜测，很难做出具有创新性的结果。开普勒的三大定律（尤其是后两条）并非一眼就能看出的，它们隐含在行星运动的轨迹之中，隐含在第谷记录下来的一大堆数据之中。历史上这样的例子实在太多了。在获得了一定数量的资料数据后，人们常常会先去猜测某些结果，然后试图去证明它。猜测一经证明就成了定理，而定理一旦插上想象的翅膀，又常常会被推广出许多更为广泛的结果。即使猜测被证明是错误的，结果也决不是一无所获的失败而常常是对问题的更为深入的了解。§2.3最短路径与最速方案问题在解决实际问题时，注意观察和善于想象是十分重要的，观察与想象例5（最短路径问题）设有一个半径为r的圆形湖，圆心为O。A、B

位于湖的两侧，AB连线过O，见图。现拟从A点步行到B点，在不得进入湖中的限制下，问怎样的路径最近。

ABOr将湖想象成凸出地面的木桩，在AB间拉一根软线，当线被拉紧时将得到最短路径。根据这样的想象，猜测可以如下得到最短路径：过A作圆的切线切圆于E，过B作圆的切线切圆于F。最短路径为由线段AE、弧EF和线段FB连接而成的连续曲线（根据对称性，AE′，弧E′F′，F′B连接而成的连续曲线也是）。EFE′F′例5（最短路径问题）设有一个半径为r的圆形湖，圆心为以上只是一种猜测，现在来证明这一猜测是正确的。为此，先介绍一下凸集与凸集的性质。定义2.1（凸集）称集合R为凸集，若x1、x2∈R及λ∈[0，1]，总有λx1+（1+λ）x2∈R。即若x1、x2∈R，则x1、x2的连线必整个地落在R中。定理2.2（分离定理）对平面中的凸集R与R外的一点K，存在直线l,l

分离R与K，即R与K分别位于l的两侧（注：对一般的凸集R与R外的一点K，则存在超平面分离R与K），见图。klR下面证明猜想以上只是一种猜测，现在来证明这一猜测是正确的。为此，先介绍一猜测证明如下：（方法一）显然，由AE、EF、FB及AE′，E′F′，F′B围成的区域R是一凸集。利用分离定理易证最短径不可能经过R外的点，若不然，设Γ为最短路径，Γ过R外的一点M，则必存在直线l分离M与R，由于路径Γ是连续曲线，由A沿Γ到M，必交l于M1，由M沿Γ到B又必交l于M2。这样，直线段M1M2的长度必小于路径M1MM2的长度，与Γ是A到B的最短路径矛盾，至此，我们已证明最短路径必在凸集R内。不妨设路径经湖的上方到达B点，则弧EF必在路径F上，又直线段AE是由A至E的最短路径，直线FB是由F到B的最短路径，猜测得证。ABOrEFE′F′M1M2MΓl猜测证明如下：（方法一）显然，由AE、EF、FB及AE′，还可用微积分方法求弧长，根据计算证明满足限止条件的其他连续曲线必具有更大的长度；此外，本猜测也可用平面几何知识加以证明等。根据猜测不难看出，例5中的条件可以大大放松，可以不必设AB过圆心，甚至可不必设湖是圆形的。例如对下图，我们可断定由A至B的最短路径必为l1与l2之一，其证明也不难类似给出。ABl1l2D到此为止，我们的研讨还只局限于平面之中，其实上述猜测可十分自然地推广到一般空间中去。1973年，J.W.Craggs证明了以上结果：若可行区域的边界是光滑曲面。则最短路径必由下列弧组成，它们或者是空间中的自然最短曲线，或者是可行区域的边界弧。而且，组成最短路径的各段弧在连接点处必定相切。还可用微积分方法求弧长，根据计算证明满足限止条件的其他连续曲例6

一辆汽车停于A处并垂直于AB方向，此汽车可转的最小圆半径为R，求不倒车而由A到B的最短路径。解（情况1）若|AB|>2R，最短路径由弧AC与切线BC