多目标规划培训教材

1、MATLABMATLAB求解多目标规划求解多目标规划江西江西师师范大范大学数学数信信学学院院吴吴根秀根秀一、一、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,

2、 exitflag = bintprog(.)x, fval, exitflag, output = bintprog(.)数学模型：数学模型：MIN lambda S.T. F(x)-weight* lambda =goal(达到目标）达到目标） Ax=b（线性不等式约束）（线性不等式约束） Aeqx=beq（线性等式约束）（线性等式约束） C(x)=0（非线性不等式约束）（非线性不等式约束） Ceq(x)=0（非线性等式约束）（非线性等式约束） lb=x 1 % Two output arguments G = . % Gradients evaluated at xEndThe grad

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

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

5、均在-4,4中中; 设特征值的设特征值的weight= abs(goal),定义目标函数定义目标函数F如下：如下： 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 matrixgoal = -5 -3 -1; % Set goal

6、values for the eigenvaluesweight = abs(goal) % Set weight for same percentagelb = -4*ones(size(K0); % Set lower bounds on the controllerub = 4*ones(size(K0); % Set upper bounds on the controlleroptions = optimset(Display,iter); % Set display parameterK,fval,attainfactor = fgoalattain(K)eigfun(K,A,B,

7、C). goal,weight,lb,ub,options)二二、举例举例-有关循环控制系统优化问题有关循环控制系统优化问题v运行结果如下运行结果如下vActive constraints:v 1v 2v 4v 9v 10vK = v -4.0000 -0.2564v -4.0000 -4.0000vfval =v -6.9313v -4.1588v -1.4099vattainfactor = v -0.3863二二、举例举例-有关循环控制系统优化问题有关循环控制系统优化问题如果至少保证如果至少保证38.63%的目标精确匹配，设置的目标精确匹配，设置GoalsExactAchieve参数值为

8、参数值为3v options = optimset(GoalsExactAchieve,3);v K,fval,attainfactor = fgoalattain(.v (K)eigfun(K,A,B,C),K0,goal,weight,lb,ub,.options)v After about seven iterations, a solution is v K = v -1.5954 1.2040v -0.4201 -2.9046v fval =v -5.0000v -3.0000v -1.0000v attainfactor = v 1.0859e-20表明目标已完全匹配表明目标已完全

9、匹配二二、举例举例-有关循环控制系统优化问题有关循环控制系统优化问题谢谢！谢谢！ 初等模型举例初等模型举例 常见类型常见类型v定性模型定性模型v经验公式（拟合、插值）经验公式（拟合、插值）v量纲分析量纲分析v比例模型比例模型 2.1 崖高的崖高的估算估算假如你站在崖顶且身上带着一只具有跑表功假如你站在崖顶且身上带着一只具有跑表功 能的计算器，你也许会出于好奇心想用扔下能的计算器，你也许会出于好奇心想用扔下 一块石头听回声的方法来估计山崖的高度，一块石头听回声的方法来估计山崖的高度， 假定你能准确地测定时间，你又怎样来推算假定你能准确地测定时间，你又怎样来推算 山崖的高度呢，请你分析一下这一问题

10、。山崖的高度呢，请你分析一下这一问题。我有一只具有跑我有一只具有跑 表功能的计算器。表功能的计算器。方法一方法一假定空气阻力不计，可以直接利用自由落体运动的公式假定空气阻力不计，可以直接利用自由落体运动的公式来计算。例如，来计算。例如， 设设t=4秒，秒，g=9.81米米/秒秒2，则可求得，则可求得h78.5米。米。221gth 我学过微积分，我可以做我学过微积分，我可以做 得更好，呵呵。得更好，呵呵。 vKmgdtdvmF除去地球吸引力外，对石块下落影响最大的当除去地球吸引力外，对石块下落影响最大的当 属属空气阻空气阻力力。根据流体力学知识，此时可设空气阻力正比于石块下。根据流体力学知识，此

11、时可设空气阻力正比于石块下落的速度，阻力系落的速度，阻力系 数数K为常数，因而，由牛顿第二定律可为常数，因而，由牛顿第二定律可得：得： kgcevkt令令k=K/m，解得解得 代入初始条件代入初始条件 v(0)=0，得，得c=g/k，故有，故有 ktekgkgv再积分一次，得：再积分一次，得： cekgtkghkt2若设若设k=0.05并仍设并仍设 t=4秒，则可求秒，则可求 得得h73.6米。米。 听到回声再按跑表，计算得到的时间中包含了听到回声再按跑表，计算得到的时间中包含了 反应时间反应时间 不妨设不妨设平均反应时间平均反应时间 为为0.1秒秒 ，假如仍，假如仍 设设t=4秒，扣除反秒，

12、扣除反应时间后应应时间后应 为为3.9秒，代入秒，代入 式式，求得，求得h69.9米。米。 222)1(kgektkgkgekgtkghktkt多测几次，取平均多测几次，取平均值值代入初始条代入初始条 件件h(0)=0，得到计算山崖高度的公式：，得到计算山崖高度的公式： 将将e-kt用泰勒公式展开并用泰勒公式展开并 令令k 0+ ，即可，即可得出前面不考虑空气阻力时的结果。得出前面不考虑空气阻力时的结果。还应考虑还应考虑回声回声传回来所需要的时间。为此，令石块下落传回来所需要的时间。为此，令石块下落 的真正时间的真正时间 为为t1，声音传回来的时间记，声音传回来的时间记 为为t2，还得解一个，

13、还得解一个方程组：方程组： 933401212211.ttthkg)ekt (kghkt这一方程组是这一方程组是非线性非线性的，求的，求解不太容易，解不太容易，为了估算崖高为了估算崖高竟要去解一个竟要去解一个非线性主程组非线性主程组似乎不合情理似乎不合情理 相对于石块速度，声音速度要快得多，我们可相对于石块速度，声音速度要快得多，我们可 用方法二先求一次用方法二先求一次 h，令，令t2=h/340，校正，校正t，求石，求石块下落时间块下落时间 t1t-t2将将t1代入式代入式再算一次，得出再算一次，得出崖高的近似值。例如，崖高的近似值。例如， 若若h=69.9米，则米，则 t20.21秒，故秒

14、，故 t13.69秒，求得秒，求得 h62.3米。米。 2.2 录像带还能录多长时间录像带还能录多长时间录像机上有一个四位计数器，一盘录像机上有一个四位计数器，一盘 180分钟分钟的录像带在开始计数时为的录像带在开始计数时为 0000，到结束时计，到结束时计数为数为1849，实际走时为，实际走时为185分分20秒。我们从秒。我们从0084观察到观察到0147共用时间共用时间3分分21秒。若录像秒。若录像机目前的计数为机目前的计数为1428，问是否还能录下一个，问是否还能录下一个 60分钟的节目？分钟的节目？rRl由由W vt)r (R22得到得到212rvtWR又又 因和因和 得得 Rl tv

15、l tRv 积分得到积分得到tdt)rWvtv(d02120rrWvtW2rWvtW2t2212021)()(即即从而有从而有rrWvtWn212)(12我们希望建立一个录像带已录像时我们希望建立一个录像带已录像时 间间t与计数器计与计数器计 数数n之间的函数关系。为建立一个正确的模型，首之间的函数关系。为建立一个正确的模型，首 先必先必须搞清哪些量是常量，哪些量是变量。首先，录像须搞清哪些量是常量，哪些量是变量。首先，录像 带带的厚的厚 度度W是常量，它被绕在一个半径是常量，它被绕在一个半径 为为r的园盘上，的园盘上，见图。磁带转动中线速见图。磁带转动中线速 度度v显然也是常数，否则图象显然

16、也是常数，否则图象声音必然会失真。此外，计数器的读声音必然会失真。此外，计数器的读 数数n与转过的圈与转过的圈数有关，从而与转过的角数有关，从而与转过的角 度度成正比。成正比。 rRlrrWvtWn212)(12 此式中的三个参数此式中的三个参数W、v和和r均不易精确测得，均不易精确测得，虽然我们可以从上式解出虽然我们可以从上式解出t与与n的函数关系，的函数关系，但效果不佳，故令但效果不佳，故令 则可将上式简化为：则可将上式简化为： Wv vW/r2tn故故nnnt21222令令21a b2上式又可化简记成上式又可化简记成 t= an2+bn t= an2+bn rRl上式以上式以a、b为参数

17、显然是一个十分明智的为参数显然是一个十分明智的做法，它为公式的最终确立即参数求解提做法，它为公式的最终确立即参数求解提供了方便。将已知条件代入，得方程组：供了方便。将已知条件代入，得方程组： 3.351471478484185.331849(1849)12122tbatbaba从后两式中消从后两式中消 去去t1，解得，解得a=0.0000291, b=0.04646,故故t=0.0000291 n2+0.04646n，令，令n=1428，得到，得到t=125.69（分）由于一盒录像带实际可录像时间为（分）由于一盒录像带实际可录像时间为185.33分，分，故尚可录像时间故尚可录像时间 为为59.

18、64分，已不能再录下一个分，已不能再录下一个60分钟分钟的节目了。的节目了。 在解决实际问题时，注意观察和善于想象是十分重要的，在解决实际问题时，注意观察和善于想象是十分重要的，观察与想象不仅能发现问题隐含的某些属性，有时还能顺观察与想象不仅能发现问题隐含的某些属性，有时还能顺理成章地找到解决实际问题的钥匙。本节的几个例子说明，理成章地找到解决实际问题的钥匙。本节的几个例子说明，猜测也是一种想象力。没有合理而又大胆的猜测，很难做猜测也是一种想象力。没有合理而又大胆的猜测，很难做出具有创新性的结果。开普勒的三大定律（尤其是后两条）出具有创新性的结果。开普勒的三大定律（尤其是后两条）并非一眼就能看

19、出的，它们隐含在行星运动的轨迹之中，并非一眼就能看出的，它们隐含在行星运动的轨迹之中，隐含在第谷记录下来的一大堆数据之中。历史上这样的例隐含在第谷记录下来的一大堆数据之中。历史上这样的例子实在太多了。在获得了一定数量的资料数据后，人们常子实在太多了。在获得了一定数量的资料数据后，人们常常会先去猜测某些结果，然后试图去证明它。猜测一经证常会先去猜测某些结果，然后试图去证明它。猜测一经证明就成了定理，而定理一旦插上想象的翅膀，又常常会被明就成了定理，而定理一旦插上想象的翅膀，又常常会被推广出许多更为广泛的结果。即使猜测被证明是错误的，推广出许多更为广泛的结果。即使猜测被证明是错误的，结果也决不是一

20、无所获的失败而常常是对问题的更为深入结果也决不是一无所获的失败而常常是对问题的更为深入的了解。的了解。 2.3最短路径与最速方案问题最短路径与最速方案问题 例例5（最短路径问题）（最短路径问题） 设有一个半径为设有一个半径为 r 的圆形湖，圆心为的圆形湖，圆心为 O。A、B 位于湖的两侧，位于湖的两侧，AB连线过连线过O，见图。，见图。现拟从现拟从A点步行到点步行到B点，在不得进入湖中的限点，在不得进入湖中的限 制下，问怎样的路径最近。制下，问怎样的路径最近。 ABOr将湖想象成凸出地面的木桩，将湖想象成凸出地面的木桩， 在在AB间拉一根软线，当间拉一根软线，当线被拉紧时将得到最短路径。根据这

21、样的想象，猜测线被拉紧时将得到最短路径。根据这样的想象，猜测 可以如下得到最短路径：可以如下得到最短路径： 过过A作圆的切线切圆于作圆的切线切圆于E，过，过B作圆的切线切圆作圆的切线切圆 于于F。最短路径为由线。最短路径为由线 段段AE、弧、弧EF和线段和线段FB连接而成的连续曲线（根据对称性，连接而成的连续曲线（根据对称性，AE，弧弧EF，FB连接而成的连续曲线也是）。连接而成的连续曲线也是）。EFEF以上只是一种猜测，现在来证明这一猜测是正确的。为此，以上只是一种猜测，现在来证明这一猜测是正确的。为此，先介绍一下凸集与凸集的性质。先介绍一下凸集与凸集的性质。定义定义2.1（凸集凸集）称集合

22、）称集合 R为凸集，若为凸集，若x1、x2R及及0，1，总有总有x1+（1+）x2R。即若。即若x1、x2R，则，则x1、x2的连线必整个地落的连线必整个地落 在在R中。中。定理定理2.2（分离定理分离定理）对平面中的凸）对平面中的凸 集集R与与R外的一点外的一点K，存在直线存在直线 l , l 分离分离R与与K，即，即R与与K分别位于分别位于 l 的两侧（注：的两侧（注：对一般的凸对一般的凸 集集R与与R外的一点外的一点K，则存在超平面分，则存在超平面分 离离R与与K），见图。），见图。klR下面证明猜想下面证明猜想猜测证明如下：猜测证明如下：（方法一）（方法一）显然，显然， 由由AE、EF

23、、FB及及AE，EF，FB围成围成的区域的区域 R是一凸集。利用是一凸集。利用分离定理分离定理易证最短径不可能经过易证最短径不可能经过R外的点，若不然，设外的点，若不然，设 为最短路径，为最短路径，过过R外的一点外的一点M，则，则必存在直必存在直 线线l分离分离M与与R，由于路径，由于路径是连续曲线，由是连续曲线，由A沿沿到到M，必交，必交l于于M1，由，由M沿沿到到B又必交又必交l于于M2。这样，直线。这样，直线 段段M1M2的长度必小于路的长度必小于路 径径M1MM2的长度，与的长度，与是是A到到B的的最短路径矛盾，至此，我们已证明最短路径必在凸集最短路径矛盾，至此，我们已证明最短路径必在

24、凸集R内。内。不妨设路径经湖的上方到达不妨设路径经湖的上方到达B点，则弧点，则弧EF必在路径必在路径F上，又上，又直线段直线段AE是由是由A至至E的最短路径，直线的最短路径，直线FB是由是由F到到B的最短的最短路径，猜测得证。路径，猜测得证。ABOrEFEFM1M2Ml还可用还可用微积分微积分方法求弧长，根据计算证方法求弧长，根据计算证明满足限止条件的其他连续曲线必具有明满足限止条件的其他连续曲线必具有更大的长度；此外，本猜测也可用更大的长度；此外，本猜测也可用平面平面几何几何知识加以证明等。知识加以证明等。 根据猜测不难看出，根据猜测不难看出， 例例5中的条件可以大大放中的条件可以大大放松，

25、可以不必松，可以不必 设设AB过圆心，甚至可不必设湖过圆心，甚至可不必设湖是圆形的。例如对是圆形的。例如对 下图，我们可断定由下图，我们可断定由A至至B的最短路径必的最短路径必 为为l1与与l2之一，其证明也不难类之一，其证明也不难类似给出。似给出。 ABl1l2D到此为止，我们的研讨还只局限于平面之中，到此为止，我们的研讨还只局限于平面之中，其实上述猜测可十分自然地推广到一般空间其实上述猜测可十分自然地推广到一般空间中去。中去。1973年，年，J.W.Craggs证明了以上结果：证明了以上结果：若可行区域的边界是光滑曲面。则最短路径必由下列弧组若可行区域的边界是光滑曲面。则最短路径必由下列弧

26、组成，它们或者是空间中的自然最短曲线，或者是可行区域成，它们或者是空间中的自然最短曲线，或者是可行区域的边界弧。而且，组成最短路径的各段弧在连接点处必定的边界弧。而且，组成最短路径的各段弧在连接点处必定相切。相切。例例6 6 一辆汽车停于一辆汽车停于 A A处并垂直于处并垂直于ABAB方向，此方向，此汽车可转的最小圆半径为汽车可转的最小圆半径为 R，求不倒车而由，求不倒车而由 A A到到B B的最短路径。的最短路径。解解（情况（情况1）若若|AB|2R，最短路径由，最短路径由 弧弧AC与切线与切线BC组组成（见成（见图图 ）。）。（情况（情况2）若若|AB|0为推力，为推力，fS，故由连续函数

27、的性质存在，故由连续函数的性质存在 某某TT，S（T）=S但这一结但这一结果与果与=(t)是最优方案下的车速的假设矛盾，因为用我们猜测是最优方案下的车速的假设矛盾，因为用我们猜测的推车方法推车，只的推车方法推车，只 需需T时间即可将车推到修车处，时间即可将车推到修车处， 而而T0可推出可推出 0的置换矩阵的置换矩阵Pijt步步2 确定确定 min|1ijijtp步步3 取取 ，用，用 代替代替PTPT步步4 若若 =0，停；否则，返回步，停；否则，返回步1。T例例2. 2. 为方便起见，我们来分解一个元素均为非负整数的为方便起见，我们来分解一个元素均为非负整数的3阶双随机矩阵，阶双随机矩阵，（

28、由（由Birkhoff定理，定理，r5）145532433T解：取解：取 ，=min 1, 3, 3 =11100010001P分解成分解成1045522432TP，再取，再取 2001100010P因因min 5, 5, 3 = 3，又有，又有120423222402TPP，取，取 3010001100P于是又有于是又有 12302232220202TPPP易得分解结果为：易得分解结果为：1000010100100010103 1002 0012 1002 010001010100001100T110rkk尚需解决的问题是如何求尚需解决的问题是如何求P，使得，使得Pij0必有必有 。读者不难

29、发现，此问。读者不难发现，此问题可以通过求解一个两分图上的最大流（或最大匹配）来实现，计算量题可以通过求解一个两分图上的最大流（或最大匹配）来实现，计算量为为O(n4)，是多项式时间可解的。具体方法为：作一两分图，若，是多项式时间可解的。具体方法为：作一两分图，若 ，则作边（则作边（i, j），令边容量为），令边容量为1，这样，可作出，这样，可作出P的充要条件是该最大流问的充要条件是该最大流问题的最大流量为题的最大流量为n。对例。对例9.33，n=3。由于所有。由于所有 ，先取，先取 0ijt 0ijt 0ijt ， P1为为 1100010001PT045522432相应的两分图为：相应的两

30、分图为：于是又可求得于是又可求得2001100010P120423222402TPP，相应的两分图为：，相应的两分图为： 又可得又可得 ，如此下去，直到作不出，如此下去，直到作不出P为至，为至，由于由于 的特殊性质及的特殊性质及Birkhoff定理，上述分解必能在不超过定理，上述分解必能在不超过r= (n1)2 + 1步内终止。步内终止。3010001100PT上述开关设计方法要求在通讯卫星上设置上述开关设计方法要求在通讯卫星上设置(n1)2 + 1种不同的开关模式种不同的开关模式（即（即Pk），当），当n稍大时，稍大时，(n1)2 + 1仍显得太大而使得使用时不便。例如，仍显得太大而使得使用

31、时不便。例如，当当n=41时，时，(n1)2 + 1=| 60 |。为实用方便，人们研究了限止开关模式个。为实用方便，人们研究了限止开关模式个数的相应问题。数的相应问题。若要求若要求rn，即要求通讯卫星上至多设置，即要求通讯卫星上至多设置n种开关模式，则问题化为令种开关模式，则问题化为令rn，求不超过求不超过n个置换矩阵个置换矩阵Pk及及k，使之满足：，使之满足： min S.t 1nkk1nkkkPT为了使任意一对发射法与接收站之间的传送均为可能实现的，自然应要求为了使任意一对发射法与接收站之间的传送均为可能实现的，自然应要求Pk满足满足（5.1） 1111111nkkP （5.2） （右面

32、的矩阵有（右面的矩阵有n2个值为个值为1的分量，每一的分量，每一Pk 恰有恰有n个个1分量）故分量）故r=n。容易看出，（容易看出，（5.1）隐含着）隐含着T的每一元素只能被唯一的的每一元素只能被唯一的P复盖，即复盖，即T的元素在的元素在分解中是不可分割的，这当然是一个好性质，使实际操作时较为方便，但分解中是不可分割的，这当然是一个好性质，使实际操作时较为方便，但可惜的是对一般的双随机矩阵，分解很可能无解。可惜的是对一般的双随机矩阵，分解很可能无解。例例3 3 若取若取1100010001P， 2010001100P3001100010P， 145532433T（注意：（注意：T已是双随机矩阵

33、，行和列和均为已是双随机矩阵，行和列和均为10） 则min 31kkS.t 31kkkPT的解为的解为1=3，2=4，3=5。3112kk（大于（大于10）而）而 31345534453kkkPT但等号经常并不成立。但等号经常并不成立。1985年，年，FRendel证明，在给定满足（证明，在给定满足（5.2）的置）的置换矩阵换矩阵P1,Pn后，求解问题（后，求解问题（5.1）是）是NP难的，从而不可能存在多项式难的，从而不可能存在多项式时间算法，除非时间算法，除非P=NP。现要求现要求r2n一种自然而方便的开关设置为引入两组各有一种自然而方便的开关设置为引入两组各有n个开关模式的置换矩阵个开关

34、模式的置换矩阵P1,Pn，Q1,Qn，满足下面的（，满足下面的（5.3)式：式：111111nnkkkkPQ例如，当例如，当n=3时，可令：时，可令：1100010001P2010001100P3001100010P1001010100Q2010100001Q3100001010Q（注：这种设置方法保持了其内在的对称性，不失为一种明智的做法。）（注：这种设置方法保持了其内在的对称性，不失为一种明智的做法。）现在，我们来分解例现在，我们来分解例9.33中的双随机矩阵中的双随机矩阵 ，令，令 =，得方程组，得方程组TT31()kkkkkPQ132231321123213312145532433求出

35、各对角线与反对角线上的三个元素之和，并作一些简单的消去运算；求出各对角线与反对角线上的三个元素之和，并作一些简单的消去运算；将矩阵的所有元素相加，可得下面的方程组：将矩阵的所有元素相加，可得下面的方程组：313212131231232102()()10注意到（注意到（5.3），易证空间），易证空间 的维数为的维数为 5，故故 之一可任取，（稍加注意即可保持非负性），之一可任取，（稍加注意即可保持非负性），例如，令例如，令3=0，求得，求得 ，故有，故有31()kkkkkPQ,kk 121232,1,2,3123122322TPPPQQ31()10kkk读者不难验证，上述方法可推广到读者不难验证

36、，上述方法可推广到n是奇数的一般情况。事实，由各对角线是奇数的一般情况。事实，由各对角线元素之和可导出元素之和可导出n1个方程，由各反对角线元素之和又可导出个方程，由各反对角线元素之和又可导出n1个方程，个方程，加上矩阵所有元素之和导出的等式，共计可导出加上矩阵所有元素之和导出的等式，共计可导出2 n1个方程，并易知它个方程，并易知它们是独立的。另一方面空间们是独立的。另一方面空间 的维数恰为的维数恰为2 n1，故，故 之一可任取，而通过方程组解得所有的之一可任取，而通过方程组解得所有的 ，（只须注意保持其非负性即可）（只须注意保持其非负性即可）1()nkkkkkPQ,kk ,kk 但当但当n

37、为偶数时，情况就不大相同了。让我们先来观察一下为偶数时，情况就不大相同了。让我们先来观察一下n=4的情况。的情况。当当n=4时，时，11000010000100001P20100001000011000P30010000110000100P40001100001000010P10001001001001000Q20010010010000001Q30100100000010010Q41000000100100100Q1423324143122134413241142321344312()kkkkKABPQBA 易见，易见， 具有非常特殊的结构，一般的偶数阶双随机矩具有非常特殊的结构，一般的偶数阶双随机矩阵，即使其元素是非负整数，也无法用阵，即使其元素是非负整数，也无法用Pk、Qk来分解。来分解。41()kkkkkPQ当当 具有上述结构时，能否用具有上述结构时，能否用Pk和和Qk来分解呢？易见，由各对角线元素来分解呢？易见，由各对角线元素之和可导出：之和可导出：T12441332421342()1042()1242()1442()16另外，由反对角线元素之和又可导出另外，由反对角线元素之和又可导出12441332421342()1442()1042()1442()14上述方程中只有上述方程中只有6个是独立的，且已不可能再得出新的独立方程，（读者个是独立的，