通过matlab的线性规划函数解决完美匹配问题

前文《GEODUAL软件说明书之完美匹配问题》中介绍了使用软件求解完美匹配问题的步骤和理论。现在通过matlab的线性规划求解器，从问题源头出发，手动求解完美匹配问题。本例采用10个点的匹配问题。坐标如图1，2所示：

完美匹配问题满足的规划问题如下：

由于护城河（moat）都是在相切圆之后确定的，因此先使护城河的宽度 为0，则上面的规划问题就变为了：

下面通过matlab来求解上述的线性规划问题。由于有10个坐标点，第一步求解点集间的距离矩阵如图5：

之后列 的约束不等式，由于有10个点，且 ，所以一共有9+8+....1=45个不等式。将其系数列入表格并且带入距离值可得：

使用matlab的linprog求解函数求解此线性规划问题的解，其中满足

求得的结果如图8所示：

此结果为每个坐标点为圆心的最大半径值，绘制出如图9所示：

根据上图可以简单地画出完美匹配问题的结果：

将其和GEODUAL软件中的结果对比，发现线段匹配的结果相同，但是相切圆的画法出现了问题，其中软件输出的B和G点相切圆半径差距较大，导致G圆无法与D圆相切从而产生 的护城河，但是根据笔者计算的结果左方的护城河集合应该是 。笔者认为，应该是线性规划所用的计算公式不同导致的，笔者的算法先省略了护城河的宽度，计算完相切圆半径后再确定护城河宽度，而GEOGUAL中是连带护城河宽度一起计算的，因此存在多组线性规划的最优解，其中一组是笔者的，G点相切圆半径大些，省去了 护城河的宽度值。而软件计算的是使G点相切圆半径小些，增加了 护城河的宽度值。而前后两种算法的增减量之和是相同的。

第一个问题比较简单，这里懒得对着你的图片敲数据，用随机数代替

long=100

A=0.1*rand(long,1)

B=0.15*rand(long,1)

[AA BB]=meshgrid(A,B)

gg=0.5017*AA-0.65*BB

g=gg>=-0.03&gg<=0.03%这就是gij矩阵，相当于二分图的连接矩阵

[num h] = maxnum(g)%第二个问题就是求二分图的最大匹配问题，这里

%调用了一个自己写maxnum函数，返回num就是最大值，h是hij(不唯一）

以下是maxnum.m的内容，用的是匈牙利算法

其中还用了一个递归的incpath函数，寻找增广路径

function [num h] = maxnum(g)

s=size(g)

global G_h%矩阵hij记录选中

global G_g%矩阵gij记录匹配

global G_v%记录当前一次路径访问过的节点

G_h=false(s)%矩阵hij初始为空

G_g=g%矩阵gij就是传递进来的参数g

for i=1:s(1)

    G_v=false(1,s(2))%每次初始化径访问过的节点为空

    incpath(i)%从Ai开始寻找增广路径

end

h=G_hnum=sum(h(:))%输出最大匹配数，和匹配矩阵h

clear global 'G_h'clear global 'G_g'

end

function OK = incpath(i)%从Ai开始

global G_hglobal G_gglobal G_vOK=false

j=find(~G_h(i,:)&G_g(i,:)&~G_v,1)%寻找合条件的Bj

if isempty(j),returnend%找不到返回false

G_v(j)=true%找到了，标记Bj为以访问节点

ii=find(G_h(:,j))%寻找Bj在原来匹配中

if isempty(ii)%如果不在原匹配中

G_h(i,j)=trueOK=truereturnend%找到增广路径末端，返回true

ok=incpath(ii)%如果在原来的匹配中，根据匹配对应的Aii递归调用incpath寻找

if ok %如果递归寻找返回成功

G_h(i,j)=~G_h(i,j)G_h(ii,j)=~G_h(ii,j)OK=truereturnend%路径反色返回true

end

我看了一下你的链接和程序.

这是你没定义dtwOptSet,当然dtw和dtwOptSet都是作者自定义的函数，不在matlab的标准库里，这个图也是明显用了3个subplot画的

如果你想运行这个，请去作者推荐的

http://mirlab.org/jang/books/dcpr/introMatlabProgram.asp?title=1-2%20Example%20Programs%20(%A6p%A6%F3%A8%FA%B1o%B5{%A6%A1%BDX)

下载example就可以了.

---