多目标规划MATLAB2012年wgx.ppt

1、MATLAB求解多目标规划,一、0-1规划的MATLAB求解,数学模型：MIN fx S.T. Ax=b Aeqx=beq x=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, outp

2、ut = bintprog(.),数学模型：MIN lambda S.T. F(x)-weight* lambda =goal(达到目标） Ax=b（线性不等式约束） Aeqx=beq（线性等式约束） C(x)=0（非线性不等式约束） Ceq(x)=0（非线性等式约束） lb=x=ub F=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

3、,x0,goal,weight,A,b,Aeq,beq) x = fgoalattain(fun,x0,goal,weight,A,b,Aeq,beq,lb,ub),二、多目标规划的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(

4、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,ex

5、itflag,output,lambda = fgoalattain(.),二、多目标规划的MATLAB求解,x = fgoalattain(myfun,x0,goal,weight,A,b,Aeq,beq,. lb,ub,mycon) where mycon is a MATLAB function such as function c,ceq = mycon(x) c = . % compute nonlinear inequalities at x. ceq = . % compute nonlinear equalities at,二、多目标规划的MATLAB求解,x = fgoala

6、ttain(myfun,x0,goal,weight) where myfun is a MATLAB function such as function F = myfun(x) F = . % Compute function values at x.,有关优化参数设置： options = optimset(GradObj,on)目标函数的梯度方向参数设置为on时，用下列函数定义： function F,G = myfun(x) F = . % Compute the function values at x if nargout 1 % Two output arguments G =

7、 . % Gradients evaluated at x End The gradient consists of the partial derivative dF/dx of each F at the point x.,二、多目标规划的MATLAB求解,二、多目标规划的MATLAB求解,有关优化参数设置： options = optimset(GradConstr,on）约束条件的梯度方向参数设置为on时，用下列函数定义： function c,ceq,GC,GCeq = mycon(x) c = . % Nonlinear inequalities at x ceq = . % No

8、nlinear equalities at x if nargout 2 % Nonlcon called with 4 outputs GC = . % Gradients of the inequalities GCeq = . % Gradients of the equalities End 注意：一般 weight=abs(goal),模型：x=(A+BKC)x+Bu，设计K满足目标： Y=Cx 1）循环系统的特征值（由命令eig(A+B*K*C)确定）的目标为goal = -5,-3,-1 2）K中元素均在-4,4中; 设特征值的weight= abs(goal),定义目标函数F如

9、下： function F = eigfun(K,A,B,C) F = sort(eig(A+B*K*C); % Evaluate objectives，由小到大排列 优化程序为： A = -0.5 0 0; 0 -2 10; 0 1 -2;B = 1 0; -2 2; 0 1; C = 1 0 0; 0 0 1; K0 = -1 -1; -1 -1; % Initialize controller matrix goal = -5 -3 -1; % Set goal values for the eigenvalues weight = abs(goal) % Set weight for

10、same percentage lb = -4*ones(size(K0); % Set lower bounds on the controller ub = 4*ones(size(K0); % Set upper bounds on the controller options = optimset(Display,iter); % Set display parameter K,fval,attainfactor = fgoalattain(K)eigfun(K,A,B,C). goal,weight,lb,ub,options),二、举例-有关循环控制系统优化问题,运行结果如下 Ac

11、tive constraints: 1 2 4 9 10 K = -4.0000 -0.2564 -4.0000 -4.0000 fval = -6.9313 -4.1588 -1.4099 attainfactor = -0.3863,二、举例-有关循环控制系统优化问题,如果至少保证38.63%的目标精确匹配，设置GoalsExactAchieve参数值为3 options = optimset(GoalsExactAchieve,3); K,fval,attainfactor = fgoalattain(. (K)eigfun(K,A,B,C),K0,goal,weight,lb,ub,.

12、options) After about seven iterations, a solution is K = -1.5954 1.2040 -0.4201 -2.9046 fval = -5.0000 -3.0000 -1.0000 attainfactor = 1.0859e-20表明目标已完全匹配,二、举例-有关循环控制系统优化问题,谢谢！,初等模型举例,常见类型,定性模型 经验公式（拟合、插值） 量纲分析 比例模型,2.1 崖高的估算,假如你站在崖顶且身上带着一只具有跑表功 能的计算器，你也许会出于好奇心想用扔下 一块石头听回声的方法来估计山崖的高度， 假定你能准确地测定时间，你又怎

13、样来推算 山崖的高度呢，请你分析一下这一问题。,方法一,我学过微积分，我可以做 得更好，呵呵。,令k=K/m，解得,代入初始条件 v(0)=0，得c=g/k，故有,再积分一次，得：,若设k=0.05并仍设 t=4秒，则可求 得h73.6米。,听到回声再按跑表，计算得到的时间中包含了 反应时间,进一步深入考虑,不妨设平均反应时间 为0.1秒 ，假如仍 设t=4秒，扣除反应时间后应 为3.9秒，代入 式，求得h69.9米。,多测几次，取平均值,再一步深入考虑,2.2 录像带还能录多长时间,录像机上有一个四位计数器，一盘 180分钟 的录像带在开始计数时为 0000，到结束时计 数为1849，实际走

14、时为185分20秒。我们从 0084观察到0147共用时间3分21秒。若录像 机目前的计数为1428，问是否还能录下一个 60分钟的节目？,又 因和 得,积分得到,即,从而有,此式中的三个参数W、v和r均不易精确测得，虽然我们可以从上式解出t与n的函数关系，但效果不佳，故令 则可将上式简化为：,故,t= an2+bn,上式以a、b为参数显然是一个十分明智的做法，它为公式的最终确立即参数求解提供了方便。将已知条件代入，得方程组：,从后两式中消 去t1，解得a=0.0000291, b=0.04646,故t=0.0000291 n2+0.04646n，令n=1428，得到t=125.69（分）由于

15、一盒录像带实际可录像时间为185.33分，故尚可录像时间 为59.64分，已不能再录下一个60分钟的节目了。,2.3最短路径与最速方案问题,例5（最短路径问题）,设有一个半径为 r 的圆形湖，圆心为 O。A、 B 位于湖的两侧，AB连线过O，见图。 现拟从A点步行到B点，在不得进入湖中的限 制下，问怎样的路径最近。,以上只是一种猜测，现在来证明这一猜测是正确的。为此，先介绍一下凸集与凸集的性质。,下面证明猜想,猜测证明如下：,还可用微积分方法求弧长，根据计算证明满足限止条件的其他连续曲线必具有更大的长度；此外，本猜测也可用平面几何知识加以证明等。,到此为止，我们的研讨还只局限于平面之中，其实上

16、述猜测可十分自然地推广到一般空间中去。1973年，J.W.Craggs证明了以上结果：,例6 一辆汽车停于 A处并垂直于AB方向，此 汽车可转的最小圆半径为 R，求不倒车而由 A到B的最短路径。,例7 驾驶一辆停于A处与AB成1角度的汽 车到B处去，已知B处要求的停车方向必须 与 AB成2角，试找出最短路径（除可转 的最小圆半径为R外，不受其他限止）。,最速方案问题,例8 将一辆急待修理的汽车由静止开始沿一 直线方向推至相隔 S米的修车处，设阻力不 计 ，推车人能使车得到的推力 f 满足: -BfA , f0为推力，f0为拉力。问怎样 推车可使车最快停于修车处。,此问题为一泛函极值问题，求解十

17、分困难，为得出一个最速方案。我们作如下猜测：,逻辑模型,例1 拟将一批尺寸为124的的商品装入尺寸为666的正方体包装箱中，问是否存在一种装法，使装入的该商品正好充满包装箱。,解 将正方体剖分成27个222的小正方体，并 按下图所示黑白相间地染色。,再将每一222的小正方体剖分成111的小正方体。 易见，27个222的正方体中，有14个是黑的，13个是白的（或13黑14白），故经两次剖分，共计有112个111的黑色小正方体和104个111的白色小正方体。 虽然包装箱的体积恰好是商品体积的27倍，但容易看到，不论将商品放置在何处，它都将占据4个黑色和4个白色的111小正方体的位置，故商品不可能充

18、满包装箱。,德国著名的艺术家Albrecht Drer(1471-1521)于1514年曾铸造了一枚名为“Melencotia I”的铜币。令人奇怪的是在这枚铜币的画面上充满了数学符号、数字及几何图形。这里，我们仅研究铜币右上角的数字问题,所谓的魔方是指由1n2这n2个正整数按一定规则排列成的一个n行n列的正方形 。n称为此魔方的阶 。,Drer魔方：4阶，每一行之和为34，每一列之和为34，对角线（或反对角线）之和是34，每个小方块中的数字之和是34，四个角上的数字加起来也是34,什么是Drer魔方,多么奇妙的魔方！,铜币铸造时间：1514年,构造魔方是一个古老的数学游戏，起初它还和神灵联系

19、在一起，带有深厚的迷信色彩。传说三千二百多年前（公元前2200年），因治水出名皇帝大禹就构造了三阶魔方（被人们称“洛书”），至今还有人把它当作符咒用于某些迷信活动，大约在十五世纪时，魔方传到了西方，著名的科尼利厄斯阿格里帕（1486-1535）先后构造出了39阶的魔方 。,如何构造魔方,奇数（不妨n=5）阶的情况,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,魔方数量随阶数n增长的速度实在是太惊人了！,同阶魔方的个数,允许构成魔方的数取任意实数,允许取实数,n阶魔方A、B，任意实数、,A+B是n阶魔方,具有指

20、定性质的魔方全体构成一个线性空间,问题已发生了实质性变化,注：刻画一个线性空间只需指出它的维数并求出此线性空间的一组基底,松驰问题的讨论,1在第一行中共有4种取法，为保持上述性质的成立，第二行中的1还有两种取法。当第二行的1也取定后，第三行与第四行的1就完全定位了，故一共可作出8个不同的最简方阵，称之为基本魔方并记之为Q1， ，Q8,仍以4阶方阵为例。 令R为行和，C为列和，D为对角线和，S为小方块和 定义0-方：R=C=D=S=0 定义1-方：R=C=D=S=4,R=C=D=S=1的方阵构成的线性空间具有什么样的性质？,类似于构造n维欧氏空间的标准基，利用0和1我们来构造一些R=C=D=S=

21、1的最简单的方阵。,显然， Drer空间（简称D空间）中任何一个元素都可以用Q1，Q2，Q8来线性表示，但它们能否构成D空间的一组基呢？,等号两边对应元素相比较，得r1=r2=r7=0， 所以 是线性无关,是D空间的最小生成集。,令D 即解方程组：,=,解得 D=,研究Albrecht Drer铸造的铜币,2004年浙江大学数学建模竞赛 （B题）通讯卫星上的开关设置 地面上存在着n个接收站与n个发送站，而在通讯卫星上则设置了若干种开关模式。开关模式可用矩阵P=(pij)来表示，若卫星可接收发送站i发射的信息并将信息传送回地面的接收站j时，矩阵中的元素pij =1，否则pij =0。通讯卫星上的

22、接收发送任务也可以用一个矩阵T=（tij）来表示，其元素tij为需经通讯卫星传递的由i发点发送到j接收点的信息量的传送时间长度。由于技术上的原因，当发送站i在发送给接收站j信息时，它不能同时发送给别的接收站信息；同样，当接收站j在接收发送站i的信息时，也不能同时接收其他发送站发送的信息。你的任务是：,设计一组开关模式，k=1, ,r（注：r应当尽可能小），使得对任意给定的任务矩阵T，卫星开关设置均能完成要求的发接收任务。 设计一个算法，在发接收任务T给出后，可根据你设计的开关模式（k=1,r）求出各开关的使用时间k，使得在完成预定传送任务的前提下使用各开关模式的总时间最短。 同样由于技术上的原

23、因，开关模式的总数r有一个上限。当需要传送的任务数数量较大时，仍无法分派任务。请你想一些办法来解决这一困难，（当然，这时你可能要作出一些牺牲，即传送时间可能会增加一些）。,问题及模型,问题的标准形式为：在地面上存在着n个收站与n个发战，而在通讯卫星上则设置了若干种开关模式。开关模式可用矩阵P=(pij)来表示，若卫星可接收发射站i发射的信息并将信息传送回地面的接收站j时，矩阵元素pij =1，否则pij =0。通讯卫星的接发任务也可用一矩阵T=（tij）来表示，其元素tij为需经通讯卫星传递的由i发点发送到j接受点的信息量的传送时间长度。问题要求求r并设计一组开关模式Pk，k=1, ,r及模式

24、Pk的使用时间k，使得在完成预定传送任务的前提下各开关模式使用的总时间最短，即要求求解下面的问题：,例1 设,这是一个有3个发送站与3个接收站的实例，tij在矩阵中已给出，例如由发站1传送到收站1的通讯量为3单位时间等。,分析 容易看出，三个发站需传送的时间分别为6、5、5；而三个收站需接收的时间分别为6、3、7。为完成全部传送任务，通讯卫星总传送时间至少应为7单位时间，即的下界为7。,由于技术上的原因，当发站i在发送给收站j信息时，它不能同时发送给别的收站信息；同样，当收站j在接收发站i的信息时，也不能同时接收其他发站发送的信息。这一要求说明，任一开关模式Pk应具有以下性质：（1）Pk的每一

25、行中有且只有一个1，每一列中也有且只有一个1；（2）所有的1均位于不同的行列中。,满足（1）、（2）的矩阵 被称为置换矩阵，n阶置换矩阵Pk共有n!个，当n较大时，我们不可能在通讯卫星上设置这么多种不同的开关模式。因而，为了设计出切实可行的开关模式，我们还得另想办法。,（设计方法1）,注意到Pk每行（或列）元素之和均为1，故不管如何指派开关的使用时间（即不论如何取k），矩阵,均具有某些特殊的性质，例如其行和（及列和）均为同一常数。这样的矩阵构成一个线性空间（参见逻辑模型第一节 Drer魔方），为减少开关模式的种类，可取此空间的一组基底作为开关模式。在使用这种开关模式时，无论T的元素tij怎么取

26、，通讯卫星对每一发（收）点的开通时间总和是恒定的。在这种开关模式下，可按如下方式指派各开关模式的使用时间：,步1 先将T改变为 ， 满足：,将T化为 的方法一般有无穷多种，如可如下化法：,令 事实上， ，（即通讯卫星传送总时间的下界）。,令,其中,用这种方法化例中的T，得到,的任一行（或列）中元素之和均为7。,定义1 称行和、列%和均相等的矩阵为双随机矩阵（Doubly stochastic matrix）,定理1 （Birkhoff定理，1944）任一n阶双随机矩阵均可写成至多 （n1）2+1个置换矩阵的非线性组合。,的分解可如下进行：,步1 选取由Pij0可推出 0的置换矩阵P,步2 确定

27、,步3 取 ，用 代替,步4 若 =0，停；否则，返回步1。,例2. 为方便起见，我们来分解一个元素均为非负整数的3阶双随机矩阵， （由Birkhoff定理，r5）,解：取 ，=min 1, 3, 3 =1,因min 5, 5, 3 = 3，又有,，取,于是又有,易得分解结果为：,尚需解决的问题是如何求P，使得Pij0必有 。读者不难发现，此问题可以通过求解一个两分图上的最大流（或最大匹配）来实现，计算量为O(n4)，是多项式时间可解的。具体方法为：作一两分图，若 ，则作边（i, j），令边容量为1，这样，可作出P的充要条件是该最大流问题的最大流量为n。对例9.33，n=3。由于所有 ，先取,

28、， P1为,相应的两分图为：,于是又可求得,，相应的两分图为：,又可得 ，如此下去，直到作不出P为至， 由于 的特殊性质及Birkhoff定理，上述分解必能在不超过r= (n1)2 + 1步内终止。,上述开关设计方法要求在通讯卫星上设置(n1)2 + 1种不同的开关模式（即Pk），当n稍大时，(n1)2 + 1仍显得太大而使得使用时不便。例如，当n=41时，(n1)2 + 1=| 60 |。为实用方便，人们研究了限止开关模式个数的相应问题。,若要求rn，即要求通讯卫星上至多设置n种开关模式，则问题化为令rn，求不超过n个置换矩阵Pk及k，使之满足：,min S.t,为了使任意一对发射法与接收站

29、之间的传送均为可能实现的，自然应要求 Pk满足,（5.1）,（5.2）,（右面的矩阵有n2个值为1的分量，每一Pk 恰有n个1分量）故r=n。,容易看出，（5.1）隐含着T的每一元素只能被唯一的P复盖，即T的元素在分解中是不可分割的，这当然是一个好性质，使实际操作时较为方便，但可惜的是对一般的双随机矩阵，分解很可能无解。,例3 若取,（注意：T已是双随机矩阵，行和列和均为10）,则min,S.t,的解为1=3，2=4，3=5。,（大于10）而,但等号经常并不成立。1985年，FRendel证明，在给定满足（5.2）的置换矩阵P1,Pn后，求解问题（5.1）是NP难的，从而不可能存在多项式时间算法，除非P=NP。,现要求r2n,一种自然而方便的开关设置为引入两组各有n个开关模式的置换矩阵P1,Pn，Q1,Qn，满足下面的（5.3)式：,例如，当n=3时，可令：,（注：这种设置方法保持了其内在的对称性，不失为一种明智的做法。）,现在，我们来分解例9.33中的双随机矩阵 ，令 = ，得方程组,求出各对角线与反对角线上的三个元素之和，并作一些简单的消去运算； 将矩阵的所有元素相加，可得下面的方程组：,注意到（5.3），易证空间 的维数为 5， 故 之一可任取，（稍加注意即可保持非负性）， 例如，令3=0，求得 ，故有,读者不难验证，上述方法可推广到n是奇数的一般情况。事实，由